0

I am trying to program in Matlab a conditional expectation of the form:

E[x|A<=x<=B] where X~N(u,s^2) (sorry, apparently the math editing here isn't what I am used to)

In Matlab, I have written up the following code:

Y=u+s*(normpdf(A,u,s)-normpdf(B,u,s))/(normcdf(B,u,s)-normcdf(A,u,s))

the problem is that it breaks down at higher values of A and B. For example, let u=0, s=1, A=10 and B=11. Simple logic says the answer should be between 10 and 11, but Matlab gives me back Inf because the denominator essentially becomes 0 while the numerator is 10^-23.

Any suggestions to make the formula give real numbers for all inputs?

Greg
  • 161
  • 6
  • Your values of u and s by definition say that 99.7% of your data lies below 3 (The 3-sigma rule). Your PDF values are therefore extremely low, and the functional form of the error function dictates the CDF value at {A,B} >> 3 will be 1. To get a decent answer for Y, A and B should be under 3*sigma. – Geodesic Jul 18 '13 at 04:20

1 Answers1

1

One way is to do the numerical integration yourself:

x = linspace(A,B,1000);
trapz(x,x.*normpdf(x,u,s)) / trapz(x,normpdf(x,u,s))

With your example values, this gives 10.0981, and it is pretty fast

Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • Numerical integration is likely to be inaccurate, given that the integrand is very nearly zero. The final result, which is the ratio of two very small, probably inaccurate, values, is likely to be inaccurate as well. Perhaps there is some analysis which could establish this approach, but on the face of it, I don't think it can be recommended. – Robert Dodier Jul 21 '13 at 19:10
  • @RobertDodier I'm not sure that the integrand being small implies inaccuracy. If that were so, both integrands could be multiplied by a large number before the computation. The graphs of the integrands actually look smooth, without any hint of numerical errors – Luis Mendo Jul 21 '13 at 19:17
  • The potential problems I see are that the value the of integrand is much less than the floating point epsilon (approximately 1e-16), and that a function (namely the exponential function) needs to be evaluated for large values -- if there are numerical accuracy problems, it is often the case that these are more pronounced the farther out in the tail one goes. That said, after evaluating the expected value with higher accuracy (bigfloat type in Maxima), I've concluded that those problems don't actually come into play. – Robert Dodier Jul 22 '13 at 01:19
  • @RobertDodier The fact that the numbers are small compared with `eps` is not a problem in itself. eps is a relative (not absolute) measurement. To see this, you can use the functional form `eps(x)` , which gives the distance from `x` to the closest number. `eps(1)` gives 2.2204e-016 (which equals `eps`), whereas `eps(1e-100)` gives `1.2690e-116`. So the distance is scaled down with x, and is always of the order of 1e-16 times x, provided x is well above `realmin`. – Luis Mendo Jul 22 '13 at 10:23