http://www.mathworks.com/matlabcentral/newsreader/view_thread/287378
MATLAB Central Newsreader  How to calculate the expectation of (UV)+
Feed for thread: How to calculate the expectation of (UV)+
enus
©19942014 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Wed, 21 Jul 2010 03:43:05 +0000
How to calculate the expectation of (UV)+
http://www.mathworks.com/matlabcentral/newsreader/view_thread/287378#764259
Hung Do
Hello,<br>
There are two independent random variables U of gamma distribution G(u,lambda1) and V of gamma distribution G(v,mu) where lambda1=0.45 and mu=1.0. Now I need to compute E[(VU)+], i.e. E[max{(VU),0}], which equals the double integral of the following function (note that mu=1.0, thus it can be omitted):<br>
<br>
fun1=@(t1,t2)(+((t2t1)>0)).*(t2t1).*(t1.^(u1)).*(t2.^(v1)).*exp((lambda1.*t1+t2)).*(lambda1.^(u))/(factorial(u1).*factorial(v1));<br>
<br>
As far as I know, there are two ways to compute this double integral in Matlab: using int and dblquad. I tried with int, but got a message saying error with "gt", so I guessed it didn't work for projection operator.<br>
Then I tried with dblquad and obtained very weird results where the expectation is decreasing in the integration ranges, for examples:<br>
u=8; v=20; <br>
dblquad(fun1,0,500,0,500)=1.3871e010<br>
dblquad(fun1,0,300,0,300)=0.7796<br>
dblquad(fun1,0,100,0,100)=1.0008<br>
and larger integration ranges resulted in "inf" errors.<br>
I would appreciate it if you guys could help me out with the following questions:<br>
1. What do you think is wrong with my caculation above? <br>
2. What should be the best way to calculate E[max{(VU),0}] in MATLAB?<br>
<br>
Thanks so much for your inputs.<br>
Hung

Wed, 21 Jul 2010 07:50:05 +0000
Re: How to calculate the expectation of (UV)+
http://www.mathworks.com/matlabcentral/newsreader/view_thread/287378#764285
Roger Stafford
"Hung Do" <dot_hung@yahoo.com> wrote in message <i25qc9$a5c$1@fred.mathworks.com>...<br>
> Hello,<br>
> There are two independent random variables U of gamma distribution G(u,lambda1) and V of gamma distribution G(v,mu) where lambda1=0.45 and mu=1.0. Now I need to compute E[(VU)+], i.e. E[max{(VU),0}], which equals the double integral of the following function (note that mu=1.0, thus it can be omitted):<br>
> <br>
> fun1=@(t1,t2)(+((t2t1)>0)).*(t2t1).*(t1.^(u1)).*(t2.^(v1)).*exp((lambda1.*t1+t2)).*(lambda1.^(u))/(factorial(u1).*factorial(v1));<br>
> <br>
> As far as I know, there are two ways to compute this double integral in Matlab: using int and dblquad. I tried with int, but got a message saying error with "gt", so I guessed it didn't work for projection operator.<br>
> Then I tried with dblquad and obtained very weird results where the expectation is decreasing in the integration ranges, for examples:<br>
> u=8; v=20; <br>
> dblquad(fun1,0,500,0,500)=1.3871e010<br>
> dblquad(fun1,0,300,0,300)=0.7796<br>
> dblquad(fun1,0,100,0,100)=1.0008<br>
> and larger integration ranges resulted in "inf" errors.<br>
> I would appreciate it if you guys could help me out with the following questions:<br>
> 1. What do you think is wrong with my caculation above? <br>
> 2. What should be the best way to calculate E[max{(VU),0}] in MATLAB?<br>
> <br>
> Thanks so much for your inputs.<br>
> Hung<br>
          <br>
You might have better luck with the quad2d function which allows either of the inner integrations limits to be functions rather than constants. For you it would be the lower limit which would be variable. This would allow you to avoid the discontinuity of the (t2t1)>0) factor. Numerical quadratures have trouble with discontinuities like this.<br>
<br>
Also you need to avoid excessively large upper limits since this can also give the quadrature process trouble. Only extend the range to include integrand values that contribute significantly to the result. Think how small the factor exp((lambda1.*t1+t2)) will be for values of t1 and t2 in the hundreds as compared with their values near zero. Matlab will think they are essentially zero "almost everywhere" and become confused. This subject has been discussed often in this newsgroup. Check out some of the pertinent threads.<br>
<br>
Roger Stafford

Wed, 21 Jul 2010 15:47:04 +0000
Re: How to calculate the expectation of (UV)+
http://www.mathworks.com/matlabcentral/newsreader/view_thread/287378#764420
Roger Stafford
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i268rd$g9j$1@fred.mathworks.com>...<br>
> You might have better luck with the quad2d function which allows either of the inner integrations limits to be functions rather than constants. For you it would be the lower limit which would be variable. This would allow you to avoid the discontinuity of the (t2t1)>0) factor. Numerical quadratures have trouble with discontinuities like this.<br>
> <br>
> Also you need to avoid excessively large upper limits since this can also give the quadrature process trouble. Only extend the range to include integrand values that contribute significantly to the result. Think how small the factor exp((lambda1.*t1+t2)) will be for values of t1 and t2 in the hundreds as compared with their values near zero. Matlab will think they are essentially zero "almost everywhere" and become confused. This subject has been discussed often in this newsgroup. Check out some of the pertinent threads.<br>
> <br>
> Roger Stafford<br>
        <br>
A further thought has occurred to me, Hung. If you make use of matlab's gammainc, the incomplete gamma function, using the 'upper' option as applied to your V distribution, you wouldn't need to do a double integration. The gammainc function will have already performed the inner integration for you. Just multiply the gammainc result using the t1 argument by the pdf of the U distribution at the same t1 and perform a single integration with respect to t1 using quad, quadl, or quadgk. Note that this also avoids the discontinuity you had earlier. Note also that quadgk supports inf as a valid limit of integration which might solve the problem you had with your upper bounds.<br>
<br>
Roger Stafford

Thu, 22 Jul 2010 05:17:04 +0000
Re: How to calculate the expectation of (UV)+
http://www.mathworks.com/matlabcentral/newsreader/view_thread/287378#764584
Roger Stafford
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i274po$3j6$1@fred.mathworks.com>...<br>
> A further thought has occurred to me, Hung. If you make use of matlab's gammainc, the incomplete gamma function, using the 'upper' option as applied to your V distribution, you wouldn't need to do a double integration. .....<br>
         <br>
In the previous article where I stated that you could use the gammainc function, I should have pointed out two things you would need to do to use it properly. First there needs to be a change of variable in using gammainc for this purpose because of the inverse scale parameters, lambda and mu. The gammainc function has no such parameters and hence the need for a change of variable. If your inner integration is taken with respect to t2, you will encounter this integral<br>
<br>
int(mu^v/gamma(v)*t2^(v1)*exp(mu*t2),t2,t1,inf)<br>
<br>
as part of your inner integral. If you substitute s2 = mu*t2, this becomes<br>
<br>
1/gamma(v)*int(s2^(v1)*exp(s2),s2,mu*t1,inf)<br>
<br>
which would be equal to gammainc(mu*t1,v,'upper'). (Of course if mu = 1, this is no change.)<br>
<br>
Second, because you are finding an expected value of VU (where it is nonnegative), you will also run into another integral on the inner part with an additional factor of t2:<br>
<br>
int(t2*mu^v/gamma(v)*t2^(v1)*exp(mu*t2),t2,t1,inf)<br>
<br>
which by the above substitution and further manipulation would turn out to be equal to<br>
<br>
v/mu*gammainc(mu*t1,v+1,'upper')<br>
<br>
The v+1 is the result of that additional t2 factor. (You had better check my manipulations here.)<br>
<br>
If on the other hand you decide to integrate with respect to t1 in the inner integral, the above considerations would apply to lambda and u, and of course you would then be using the 'lower' option in gammainc.<br>
<br>
As I said earlier, all this has to do with using gammainc in place of an inner integration. You will still have to call upon a numerical quadrature routine, preferably quadgk, to numerically integrate the resulting function you obtain from the above as the outer integral.<br>
<br>
Roger Stafford