Each of the two integrals (g1, g2) seem to be divergent (or at least is
considered to be so by R :) )
Try this:
z <- c(80, 20, 40, 30)
"f1" <- function(x, y, z) {dgamma(cumsum(z)[-length(z)], shape=x, rate=y)}
"g1" <- function(y, z) {integrate(function(x) {f1(x=x, y=y, z=z)}, 0.1,
0.5, rel.tol=1.5e-20)$value}
"g2" <- function(x, z) {integrate(function(y) {f1(x=x, y=y, z=z)}, 0.1,
0.5,rel.tol=1.5e-20)$value}
integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
Error in integrate(function(x) { : the integral is probably divergent
integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
Error in integrate(function(x) { : the integral is probably divergent
Are you sure:
a) you are not integrating over a singularity (in that case, you'd have to
calculate the integral as a Cauchy principal value)
b) you are not over/underflowing during your "manual integration" of the Gamma
density (would R's integrate function report this? - I don't know)
c) there is a loss of precision during numerical integration inside the g1/g2
functions
To see that c) is indeed possible try your code with z set to
1) c(8,2,4,3)*0.5 :
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 0.002982671
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 0.002985362
2) c(8,2,4,3) :
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 0.0006453548
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 0.0006466886
3) c(8,2,4,3)*5:
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 1.198051e-06
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 1.253991e-06
4) c(8,2,4,3)*20:
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 5.832236e-13
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 3.230487e-13
to get an idea of how this would play out.
You could try to tweak integrate (e.g. increasing the maximum number of
subdivisions, changing the absolute / relative precision) but
your integrand (? are you trying to marginalize over nuisance parameters in a
Bayesian setting) is numerically ill-behaved that I do not think
you may be able to integrate it in a straightforward manner.
You may try want to try:
a) adapt (from library adapt) to carry out the integration over the shape and
scale simultaneously
b) evaluate one or (even both!) integrations symbolically if possible (by hand,
tables, or CAS e.g. Mathematica or Maxima) and recode accordingly
Click on the following link to symbolically evaluate the indefinite numerical
integral of the gamma density w.r.t to its scale parameter:
http://integrals.wolfram.com/index.jsp?expr=y^(a-1)*Exp[-y%2Fx]%2F(x^s*Gamma[a])&random=false
This will give you an idea of the numerical nastiness that the poor integrate
function has to deal with .. :)
_________________________________________________________________
Show them the way! Add maps and directions to your party invites.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
