Integrating small valued functions

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!

[Torsten]

Maybe I'm wrong, but in my opinion, the expected biggest number among N random numbers generated by a normal (m,s) is given by

integral_{-Inf}^{+Inf} x*n*cdfnorm(x,m,s)^(n-1)*pdfnorm(x,m,s) dx

and can be evaluated without problems (at least for values of n up to 1000).

Best wishes

Torsten.

[Torsten]

I just noticed that our expressions for the expected biggest number are identical since

integral_{-Inf}^{+Inf} n*cdfnorm(x,m,s)^(n-1)*pdfnorm(x,m,s) dx = 1

Best wishes

Torsten.

[John D'Errico]

I predict you will find this almost impossible to get useful results from this for large values of n, at least in double precision arithmetic. You can try, but expect it to fail. Instead, you will need to work using tools that can support higher precision, or do what you can with the symbolic TB.

[Torsten]

I evaluated your integral. It is given by

m-1/Sqrt[Pi]*n*Sqrt[2]*s*integral_{-Inf}^{Inf} x*(0.5*Erfc[x])^(n-1)*Exp[-x^2] dx

and the last expression can be integrated up to high values of n.

Best wishes

Torsten.

[Paulo F.]

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

Answer 1

[Mike Hosea]

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.]

Thanks for your proposal, works beautifully.