Prev: Removing scaling label from axis?
Next: MuPAD: straightforward simplification
From: Hung Do on 20 Jul 2010 23:43
Hello,
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[(V-U)+], i.e. E[max{(V-U),0}], which equals the double integral of the following function (note that mu=1.0, thus it can be omitted):

fun1=@(t1,t2)(+((t2-t1)>0)).*(t2-t1).*(t1.^(u-1)).*(t2.^(v-1)).*exp(-(lambda1.*t1+t2)).*(lambda1.^(u))/(factorial(u-1).*factorial(v-1));

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.
Then I tried with dblquad and obtained very weird results where the expectation is decreasing in the integration ranges, for examples:
u=8; v=20;
dblquad(fun1,0,500,0,500)=1.3871e-010
dblquad(fun1,0,300,0,300)=0.7796
dblquad(fun1,0,100,0,100)=1.0008
and larger integration ranges resulted in "inf" errors.
I would appreciate it if you guys could help me out with the following questions:
1. What do you think is wrong with my caculation above?
2. What should be the best way to calculate E[max{(V-U),0}] in MATLAB?

Thanks so much for your inputs.
Hung
From: Roger Stafford on 21 Jul 2010 03:50
"Hung Do" <dot_hung(a)yahoo.com> wrote in message <i25qc9$a5c$1(a)fred.mathworks.com>...
> Hello,
> 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[(V-U)+], i.e. E[max{(V-U),0}], which equals the double integral of the following function (note that mu=1.0, thus it can be omitted):
>
> fun1=@(t1,t2)(+((t2-t1)>0)).*(t2-t1).*(t1.^(u-1)).*(t2.^(v-1)).*exp(-(lambda1.*t1+t2)).*(lambda1.^(u))/(factorial(u-1).*factorial(v-1));
>
> 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.
> Then I tried with dblquad and obtained very weird results where the expectation is decreasing in the integration ranges, for examples:
> u=8; v=20;
> dblquad(fun1,0,500,0,500)=1.3871e-010
> dblquad(fun1,0,300,0,300)=0.7796
> dblquad(fun1,0,100,0,100)=1.0008
> and larger integration ranges resulted in "inf" errors.
> I would appreciate it if you guys could help me out with the following questions:
> 1. What do you think is wrong with my caculation above?
> 2. What should be the best way to calculate E[max{(V-U),0}] in MATLAB?
>
> Thanks so much for your inputs.
> Hung
- - - - - - - - - - -
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 (t2-t1)>0) factor. Numerical quadratures have trouble with discontinuities like this.

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.

Roger Stafford
From: Roger Stafford on 21 Jul 2010 11:47
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <i268rd$g9j$1(a)fred.mathworks.com>...
> 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 (t2-t1)>0) factor. Numerical quadratures have trouble with discontinuities like this.
>
> 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.
>
> Roger Stafford
- - - - - - - - -
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.

Roger Stafford
From: Roger Stafford on 22 Jul 2010 01:17
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <i274po$3j6$1(a)fred.mathworks.com>...
> 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. .....
- - - - - - - - -
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

int(mu^v/gamma(v)*t2^(v-1)*exp(-mu*t2),t2,t1,inf)

as part of your inner integral. If you substitute s2 = mu*t2, this becomes

1/gamma(v)*int(s2^(v-1)*exp(-s2),s2,mu*t1,inf)

which would be equal to gammainc(mu*t1,v,'upper'). (Of course if mu = 1, this is no change.)

Second, because you are finding an expected value of V-U (where it is non-negative), you will also run into another integral on the inner part with an additional factor of t2:

int(t2*mu^v/gamma(v)*t2^(v-1)*exp(-mu*t2),t2,t1,inf)

which by the above substitution and further manipulation would turn out to be equal to

v/mu*gammainc(mu*t1,v+1,'upper')

The v+1 is the result of that additional t2 factor. (You had better check my manipulations here.)

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.

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.

Roger Stafford
 | 
Pages: 1
Prev: Removing scaling label from axis?
Next: MuPAD: straightforward simplification