Numerical Expectation of a Compiled function

Question

I have an integrand built as a compiled function. A simple test function is:

  integrand = Compile[{{x, _Real, 0}, {p, _Real, 0}}, 
           x^2 + p,
           RuntimeOptions -> {"EvaluateSymbolically" -> False} ];

I want to compute the expected value of this integrand using the function NExpectation, with

 NExpectation[integrand[x, 1], 
     x \[Distributed] NormalDistribution[5, 1], 
     Method -> {Automatic, "SymbolicProcessing" -> 0}]

The reason I am not using NIntegrate directly is that I'm hoping to obtain some gains in the timing (this integration will be repeated a large number of times) since I know that the density will always be a Normal.

I always get a number as an output; however, I always get before a bunch of

 CompiledFunction::cfsa: "Argument 5+x at position 1 should be a \!  
  (\"machine-size real number\"\)."

The "Symbolic" options are an attempt at eliminating the problem, which obviously don't work. I would like to eliminate the source of the error, since this is slowing down the whole computation.

[Edit:] A suggestion in the comments has been to use Monte Carlo methods for integration, as:

 NExpectation[integrand[x, 1], 
     x \[Distributed] NormalDistribution[5, 1], 
     Method -> {"MonteCarlo"}]

This solution takes care of the problem unless I want to control the number of points. This control is important since, in tests I have run, the third decimal of the integral already differs from the same integral computed in slower but more precise ways. By running

 NExpectation[integrand[x, 1], 
     x \[Distributed] NormalDistribution[5, 1], 
     Method -> {"MonteCarloRule","Points"->1000}]

the error message returns. Are there any other options?

Thanks.

Answer 1

Below is a test of the hypotheses in the comments:

The integrand function consists of the product of several Normal CDFs (different from each other). As the problem scales up, we'd have 1) a few thousand of these Normal CDF in the integrand, and 2) an equal number of integrations to be computed. -- Fred Dec 8 at 4:23

It might be faster to use an integration rule directly. – Michael E2 Dec 7 at 20:40

This is a proof-of-concept answer. Without actual functions to work on, the problems of (1) limiting the domain of integration and (2) determining the order of the integration rule remain to be worked out. It should be possible to solve these problems for cases like the following example. While it's fairly simple to solve them in this case, a robust solution would have a fast way to set up the integral in all cases. In the OP's case, situation may lie somewhere between, and further information would be welcome.

Features:

• The basic integration can be done in about 0.0025 sec. if the integration nodes and weights have been precomputed, which should be possible. If not, it takes about 0.004 sec., assuming the number of nodes needed stays roughly constant, around 50-60.
• Solving for the end points takes about 0.12 sec. It is probably possible to do this once per set of CDFs.
• The solution is about as close to machine precision accuracy as I could get (~14-15 digits). A lower precision goal should make things easier.

Example

Set-up: 2000 normal CDFs with means and standard deviations, each uniformly distributed in an interval.

pdf0 = PDF[NormalDistribution[m, s], x];   (* symbolic pdf *)
cdf0 = CDF[NormalDistribution[m, s], x];   (* symbolic cdf *)
dim = {2000};

SeedRandom[0]; (* for reproducibility *)
means = RandomReal[{4, 6}, dim];  (* means and s.d.s for CDFs *)
sds = RandomReal[{2, 5}, dim];
m0 = 5;                           (* mean and s.d. for PDF *)
s0 = 1;

Integrand: The product of the normal PDF and the CDFs.

Block[{x, m0, s0, m, s},
  With[{pdf = (pdf0 /. {m -> m0, s -> s0}), cdf = cdf0},
   integrand = 
    Compile[{{x, _Real}, {m0, _Real}, {s0, _Real}, {m, _Real, 1}, {s, _Real, 1}},
     pdf*Apply[Times, cdf],
     RuntimeAttributes -> {Listable}, Parallelization -> True
     ]
   ]];

Integration: For the purpose of analysis, the "exact" value was computed with a high-order Gauss-Kronrod rule over a larger interval that had an error estimate of 3.6*10^-43. This need not be done (and the code is not shown).

{abscissas, weights, errorweights} = NIntegrate`GaussRuleData[60, MachinePrecision];
Block[{x1 = Max[means, m0] + Max[sds, s0], x2 = Max[means, m0] + 13.997 s0},
  abscissas = Rescale[abscissas, {0, 1}, {x1, x2}];
  weights = (x2 - x1)*weights;
  weights.integrand[abscissas, m0, s0, means, sds]
  ] // AbsoluteTiming
Last@% - exact
(*
  {0.002433, 1.79187*10^-27}
  7.17465*10^-43^-43
*)

Notes & Analysis

First, let's look at what we're integrating. Below is a plot of the exponent (base 10) of the integrand and pdf. The horizontal part in blue is where there is underflow in integrand at the level of Log10[$MinMachineNumber] == -307.653.

Fig. 1. The plot of the exponent (base 10) of the integrand (blue) and the PDF (yellow).

Plot[{RealExponent@integrand[x, m0, s0, means, sds], 
  RealExponent@PDF[NormalDistribution[m0, s0], x]}, {x, -10, 50}]

The integrand is a highly rapidly vanishing, analytic function. It is possible to limit the interval of integration to a finite interval and suffer no loss of precision, when computing the integral at a given finite precision (such as machine precision above). Clearly we can limit the interval to where there is no underflow, but we can do better. Once the value of the integrand falls below its maximum by a couple orders of magnitude, the error in truncating the integration at that point will be less than the value of the function.

To figure out such points, one can use functions like FindMaximum and FindRoot. Note that the product of the CDFs is nearly equal to zero on left and nearly equal to one on the right. So the integrand is given by the PDF sufficiently far to the left. What happens on the right (that is significant to the value of the integral) depends on the relation of the mean of PDF to the means of the CDFs. If the mean of the PDF is sufficiently greater than those of the CDF, the PDF will dominate.

Side note on the limits of integration: The limits shown above were found in a cheating sort of way. As mentioned in the introduction, they can be found in a more automatic way. The lower limit can be found as follows:

(* test starting points *)
integrand[Max[means, m0], m0, s0, means, sds]
integrand[Max[means, m0] + Max[sds, s0], m0, s0, means, sds]
(*
  0.
  2.51444*10^-54
*)

The upper limit was found through a quickly estimate of the integral 1.624`*^-27and then with the precision factor 1.*^-17, the following estimates the coefficient of the s.d. s0:

Sqrt[2] InverseErfc[1.624`*^-27 * 1.*^-17]
(*  13.9971  *)

Summary

If the same set of CDFs are to be re-used a couple thousand times, I think all the integrals could be set up with a one-time code that would take much less than a second. Each integral would then take only a few milliseconds.