I have this joint beta distribution I am trying to integrate over:
alpha1 = 400
beta1 = 26000
alpha2 = 410
beta2 = 26000
integrate(function(y) {
sapply(y, function(y) {integrate(function(x) dbeta(x, shape1 = alpha1,
shape2 = beta1)*dbeta(y, shape1 = alpha2, shape2 = beta2), y,1)$value})
}, 0, 1)
This works for smaller alpha and beta even up to these values. But if I go much larger in beta the integrate function starts breaking. I've been comparing this to the Monte Carlo integration
n_sim = 1000000 # number of simulations
y = rbeta(n_sim, shape1 = alpha2, shape2 = beta2)
C2 = (1-pbeta(y, shape1 = alpha1, shape2 = beta1))
mean(C2)
which yields approximately the same answers for smaller reasonable alpha and beta.
And yes I need larger beta, I understand questioning why I need such large beta. It follows from large data sets. Does anybody know why this breaks? Maybe the inner workings of the integrate function? Is there a work around?