topic Re: error in evaluating a double integral in SAS/IML Software and Matrix Computations

error in evaluating a double integral

Emara — Fri, 26 Jul 2013 15:10:30 GMT

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

Re: error in evaluating a double integral

Rick_SAS — Fri, 26 Jul 2013 15:32:20 GMT

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 */

Re: error in evaluating a double integral

Emara — Sat, 27 Jul 2013 11:45:34 GMT

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.

Re: error in evaluating a double integral

Emara — Sun, 28 Jul 2013 21:42:41 GMT

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

Re: error in evaluating a double integral

Rick_SAS — Sun, 28 Jul 2013 23:16:49 GMT

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

Re: error in evaluating a double integral

Emara — Mon, 29 Jul 2013 14:01:53 GMT

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

Re: error in evaluating a double integral

SteveDenham — Mon, 29 Jul 2013 18:57:00 GMT

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

Re: error in evaluating a double integral

Rick_SAS — Mon, 29 Jul 2013 19:42:17 GMT

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.

Re: error in evaluating a double integral

SteveDenham — Mon, 29 Jul 2013 20:00:00 GMT

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

Re: error in evaluating a double integral

Rick_SAS — Tue, 30 Jul 2013 18:46:24 GMT

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.

Re: error in evaluating a double integral

Emara — Tue, 30 Jul 2013 20:50:03 GMT

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

Re: error in evaluating a double integral

Rick_SAS — Tue, 30 Jul 2013 21:07:08 GMT

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.

Re: error in evaluating a double integral

Emara — Tue, 30 Jul 2013 22:12:52 GMT

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

Re: error in evaluating a double integral

Rick_SAS — Fri, 09 Aug 2013 12:47:20 GMT

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.

Re: error in evaluating a double integral

Emara — Sun, 18 Aug 2013 12:50:35 GMT

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

Re: error in evaluating a double integral

Rick_SAS — Sun, 18 Aug 2013 20:06:08 GMT

I'm sorry for the delay, but I have a lot of travel in August.

I don't know whether it will work, but I hoped you might be able to scale the integral as follows.

Let h be the integrand. Then \integral h(x) dx = h(0) * \integral (h(x)/h(0)) dx.

Let g(x) = h(x) / h(0). My hope was that the integral of g is better behaved an remains bounded as your parameter varies. The value of the integral should be small for all values of the parameter. However, I haven't had time to think about this problem recently or to see whether my idea actually works.

A double integral of a convolution is just too hard for a casual observer to offer much specific help. It would take me quite a bit of effort to understand the details of your problem. I still think your best bet is still to simplify the problem and understand each component. The key is understanding why the integral isn't converging. For example, I notice that the documentation example for the QUAD function says that the inner integral should not range over an infinite domain. Have you tried integrating the inner intergral on a finiate domain such as [-10, 10]?

Re: error in evaluating a double integral

Emara — Mon, 19 Aug 2013 20:18:09 GMT

Sorry Rrick for troubling you, I really appreciate your continious efforts. I tried your first suggestion "scaling the integral". As I understood, I should substitute with ZED=0 and QQ=0 in the integrand. However, this substitution yields h(0)=0. Am I correct?

Aya

Re: error in evaluating a double integral

Emara — Mon, 19 Aug 2013 20:41:54 GMT

Hello Rick,I tried the following bounds for the inner integral [-10,10] . I got errors concerning matrix inverse and convergence. The log file is attached