How to solve this definite integral in MATLAB?!!

Hi there,

I want to solve the equation attached but something's going wrong! Please note that the Capital phi and small phi stand for "CDF" and "PDF" of a standard normal variate, respectively. I have written the following but it doesn't work! Could any of you help me please?

A = 1 - int(((normcdf((3.0+sqrt(0.85)*x)/(sqrt(1-0.85))))^7)*(normpdf(x)), x=-inf..inf)

Answer 1

[Walter Roberson]

There is no closed form solution.

If you replace the 3 by A, and the 0.85 by B (to get the general form), then about the most compact form of the answer is, in MuPAD format,

1-(int((1/256)*(1+erf((1/2)*(A+B^(1/2)*x)*2^(1/2)/(1-B)^(1/2)))^7*2^(1/2)*exp(-(1/2)*x^2)/PI^(1/2), x = -infinity .. infinity))

To get a numeric result for this with A = 3 and B = 0.85, you need to do numeric integration. The result is approximately 0.0047768095540396

Answer 2

[John BG]

Hi Mamali

what is one of the most useful techniques of probability?

if you don't know, assume uniform distribution

Which is it all about probability, how to walk through a dark room without knowing where the furniture is, if you understand what I mean.

So, you brought a CDF and pdf without any further information: again, let if be uniform.

pdf for uniform is as flat as it gets, and CDF of a flat constant is? yes, t.

had a first look to this monster function:

t=[-100:.01:100]
k1=sqrt(.85);k2=sqrt(1-.85);k3=3;x=((k3+k1*t)/k2).^7
plot(t,x)
grid on
grid minor

both sides diverge but it looks more or less symmetric, so it may be integrable. replace pdf(t)=1 and let me get to the CDF a bit later.

What would happen if you integrate the function without the CDF (big PHI)?

    y=@(t) ((k3+k1*t)/k2).^7
    y = 
        @(t)((k3+k1*t)/k2).^7
    >> integral(y,-Inf,Inf)
    Warning: Minimum step size reached near x = -Inf. There may be
    a singularity, or the tolerances may be too tight for this
    problem. 
    > In integralCalc/checkSpacing (line 456)
      In integralCalc/iterateScalarValued (line 319)
      In integralCalc/vadapt (line 132)
      In integralCalc (line 103)
      In integral (line 88) 
    ans =
     13226991826182229000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.00

to get more precision, need the Symbolic tool box, and use function vpa, variable precision arithmetic.

However, MATLAB already warning that there may be a singularity, and if not knowing anything else of that odd spots where integral may walk through, or may stuck before being able to carry on integrating towards Inf, there isn't much more that can be done.

Now the CDF of a uniform distribution, that has the pdf multiplying the function to integrate. uniform[0 1] or [-A A] makes sense, because how do you apply the concept of probability to an infinite set? you need to define a subset to then define the probabilities.

uniform pdf is a pulse, from any 2 fixed values[a b]. The value of the uniform pdf is 1/abs(a-b) throughout the interval where the pdf is defined, AND ZERO ELSEWHERE. For convenience I define the interval [0 1] pdf is 1 between 0 and 1.

The CDF of a uniform pdf [a b] is a steady ramp from zero on a to 1 on b. Here the integral interval has been reduced to [0 1] therefore I would rid the CDF (big PHI) by replacing it with a simple multiplying t.

Through the above reasonable simplifications, your reference vector t, your function x, and the integral, now are:

t=[0:.001:1]
x=@(t) t.*((k3+k1*t)/k2).^7
integral(x,0,1)
3307946.04

good to go?

if they keep throwing at you these kind of functions, it would be greatly appreciated the distributions to be defined in advance. Hope it helps John [e-mail]

[Walter Roberson]

The poster did say "standard normal variate", so assumption of uniform is mistaken.

[John BG]

you are right the question mentions standard normal variate, but the 'normality' comes from the normcdf used to jam the integral.

Assuming normal distribution of the variable inside the normcdf and a bit more allows to turn the normcdf() into a ramp-like multiplier that in turn allows to solve the integral.

If you have an alternative way, please go ahead and show us.

John

[Walter Roberson]

I already posted the formula based upon the user's own use of normcdf() and normpdf().

In order to get such a large positive output from the formula that has 1 minus the integral, your integral must have been negative. However, if you look at the formula for normpdf() you can see that it is obviously never negative (unless the standard deviation is negative!), and normcdf() is a probability and so is never negative. The 7th power of a nonnegative value is non-negative, and multiplied by a non-negative value is going to have to have a non-negative result. Therefore the integral is of values that are everywhere non-negative, so the result of the integral must be non-negative. This analysis is not enough to determine whether the integral could ever be greater than 1 (which would give a negative overall result) but it does establish that your calculation must be incorrect. (Or did you not take the "1-" into account?)

[John BG]

the 1-? no i did not miss it, i focussed on the integral and left the easy part for you.

John