Integrating small valued functions

1 次查看(过去 30 天)

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!

  5 个评论
Torsten
Torsten 2015-2-20
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.
Paulo F. 2015-2-20
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
Mike Hosea 2015-2-20
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;

更多回答(0 个)

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by