error in evaluating a double integral

[Emara]

Dear all,

I am evaluating a double integral using the quad function but i got some errors. The integral is f(Z).g(Q).T(Z,Q) where Z is N(0,1) and Q follows a scaled chi distribution .The program is attached and here is a part of the error:" Convergence could not be attained over the subinterval ( .M , .P)". Could you help me to figure out the problem?

Thanks

Aya

[Rick_SAS]

Be sure to test and debug each function (integrand) before you try to integrate them.  Graphing the functions are a good way to ensure that they are implemented correctly.

An integrand will return a scalar quantity when it receives a scalar input.  It looks like the INNER function is returning a vector of values:

Q = 1;

y = inner(0);

print y; /* should be a scalar */

[Emara]

Thanks alot Rick. The inner function was returning a vector of values. I corrected this mistake, However I got other errors. I don't know how could I overcome these errors. I think this time it is related to the singularity of a matrix!! The modified program is attached.

[Emara]

I am sorry that I am asking too many questions concerning this code.But it is really very important in my research.   I have been trying to correct the errors concerning the singularity of a matrix, so i used ginv function to avoid these errors. However, I got other errors related to

Convergence could not be attained over the subinterval

             ( 0 , .P)

Convergence could not be attained over the subinterval

             ( .M , .P)

Could you please help me to point out the reasons behind these errors?

Waiting for your reply,

Aya     

[Rick_SAS]

Have you already graphed the INNER function for various parameter values and determined that it is correct?

[Emara]

The plots are attached. Each plot is for a certain QQ value (1,3,5) and different ZED values(-5,-4,-3,-2,-1,0,1,2,3,4,5). I have made many trials to have this integration done, but I am always faced with errors. Sometimes related to the singularity of a matrix, so i tried the " ginv"  function, so I got errors related to the convergence of the integral over the specified interval. I really can't figure out what the problem is. I tried running the syntax for other integration limits; ZEDinterval (-1,1) and Qinterval (0,1), It gave me an answer  without errors " 1.659E-59". The correct answer over the whole interval (-inf,inf) and(0,inf) is supposed to be 339.4!!

Aya

[SteveDenham]

I think the convergence problem may be due to scaling.  The y axis ranges from near 0 to 10^15.  I would suggest integrating the log of the function, and applying the chain rule to get back to the integrated value of the basis function.

Steve Denham

[Rick_SAS]

OK, Steve, I'll bite. I don't see how you propose to do this. How does finding Integral Log(f(x)) dx give us Integral f(x) dx?  Furthermore, the OP hasn't claimed that f(x) > 0, which makes log(f(x)) problematic.

[SteveDenham]

While the OP doesn't say f(x) > 0, the plots that were attached certainly were defined only on f(x) > 0.  So, here's my thinking.  We want Integral f(x) dx.  We find Integral log(f(x)) dx, by taking the log of the response variable, and integrating.  Then I was thinking because Integral (log x) dx was x log x - x, we could just substitute.  But after putting pencil to paper, that's probably not going to work, as the the original function won't be unaffected by the integration.

Steve Denham

[Rick_SAS]

It's always hard to understand complicated code that someone else wrote, but at the end of the INNER module you have

   normal=(1/sqrt(2*pi))*exp(-0.5*(ZED**2));

   b=(v*avrl)cc*normal;

   return(b);

To me, it looks like NORMAL is simply the standard normal PDF.  The quantity c=(v*avl) is a constant that depends on Q.  Therefore the

integral of INNER over the domain (-Infinity, Infinity) is c, because the integral of the PDF is unity.  Therefore it is not necessary to call the QUAD function for INNER.

It looks like OUTER might be handled similarly.  If OUTER is some multiple k of the chi-square distribution, then the integral on [0, infinity) is simply k.

Incidentally, if you are interested in other limits of integration, you can use the CDF function to compute integrals of a PDF.

[Emara]

Thank you so much for your efforts. i really appreciate it.  My integration is the product of three quantities

1-standard normal pdf "normal"

2- a scaled chi distribution with degrees of freedom "df"

3- The quantity (v*avrl) which is a function of both variables QQ and ZED where ZED is the variable following the standard normal distribution and QQ is the variable following a scaled chi distribution.

The inner module is used to compute the quantity "v*avrl".  After computing the quantitty "v*avrl" which depends on both ZED and QQ, the inner integrand is computed by multipying "v*avrl" by the normal density.

The result of the inner integration should be a function of QQ which is called "w" in our syntax. This w is multiplied by the scaled chi distribution: (1/cf)*(1/sqrt(df))*(2**(1-(df/2)))*(QQ**(df-1))*(exp((-QQ**2)/2))/gamma(df/2);

Thus, i think I have to call QUAD for both the inner and outer integrals since "v*avrl" depends on both QQ and ZED, see this part in the syntax:

QQstar1=(QQ*z[i, ])+(ZED/(sqrt(mm)));

QQstar=QQstar1-shift;

kone=QQstar+QQ*fun1;

ktwo=QQstar+QQ*fun2;

pij=probnorm(kone)-probnorm(ktwo);

Thus "v*avrl"is not a constant with respect to the inner integral. Am I correct??

Aya

[Rick_SAS]

Yep. I missed that v*avrl depends on ZED. Too bad.

I don't suppose you know anything analytically about the function f(z)=v(z)*avrl(z), such as  values at +/- infinity and derivative info? If you do, you might be able to integrate by parts using the fact that phi(x) has a known antiderivative, the normal CDF.

[Emara]

No, I don't know the antiderivative information. What I know is that the v*avrl is the average run length of a certain control chart. I am computing the average run length using the markov chain technique. The martrix "r" is the transitional probability matrix and I have t  transiet states.

Aya

[Rick_SAS]

I recommend the following: factor out  c = f(0)=v(0)*avrl(0) from the integral.  Then you have the integral of

f(z)*phi(z) where phi is the standard normal PDF.  By scaling the integral, you might get convergence.

You can call QUAD on this bounded integrand, and then multiply the result of QUAD by c=f(0) in the optimization code.

[Emara]

Sorry for repeating my question, but I don't know what do you mean by scaling the integral? I tried to understand the procedure you explained but I couldn't get the point.

Thanks,

Aya Emara