# Integrating small valued functions

5 views (last 30 days)
Paulo F. on 17 Feb 2015
Commented: Paulo F. on 20 Feb 2015

Well, maybe it is better to start with the mathematical part of the problem:

I want to obtain the expected highest number among N random numbers generated by a Normal (m,s).

The expression for this would be (in LaTeX notation), f() and F() are density and cumulative functions respectively:

$$\frac{\int_{-\infty}^\infty x f(x|m,s) F^{n-1}(x|m,s)dx}{\int_{-\infty}^\infty f(x|m,s) F^{n-1}(x|m,s)dx}$$

Now, for big moments (m,s) this gives crazy values, and big moments I refer to things like 30 and above. Also crazy things start to happen with values of N bigger than 35.

An example of what happens is giving a result below the mean, when it is easy to proof that this number will always be bigger than the mean (note F()<1, => F()^k = over-representing the numbers to right of the distribution).

In code, what I do is to set up the parameters and then define the function:

f1=@(x)x.*pdfnorm(x,m,s).*cdfnorm(x,m,s).^(n-1);
f1s=@(x)pdfnorm(x,m,s).*cdfnorm(x,m,s).^(n-1);
p1=integral(@(x)f1s,-Inf,Inf);
p2=integral(@(x)f1,-Inf,Inf);
result=p2/p1;


I have read that Integral might have some troubles with very small or big numbers, as it would be the case for $x$ at some point, or cdfnorm()^{n-1} with high $n$ values. I have not find a way to solve this, do you have any clue?.

I thank you all in advance,

Cheers!

Paulo F. on 20 Feb 2015
Thanks for your answers, I will try to do it in the expression instead of using the pdf and cdf functions to see if it works. I will let you know.
Thanks again

Mike Hosea on 20 Feb 2015
Could you try this and tell me what values of m, s, and n lead to unexpected results? I added some waypoints near the center of the distribution for fear that sometimes the integrator might "miss the action" with its initial mesh.
function result = testfun(m,s,n)
f1=@(x)x.*normpdf(x,m,s).*normcdf(x,m,s).^(n-1);
f1s=@(x)normpdf(x,m,s).*normcdf(x,m,s).^(n-1);
w = m + s*linspace(-3,3);
p1=integral(f1s,-inf,inf,'Waypoints',w);
p2=integral(f1,-inf,inf,'Waypoints',w);
result=p2/p1;
Paulo F. on 20 Feb 2015
Thanks for your proposal, works beautifully.