Hello,

We are trying to use R to simulate a model based on parameters 'a' and 'b'.
This involves the following integration:

model<-function(s,x,a,b)(exp(-s*x*10^-5.5)*(s^(a-1)*(1-s)^(b-1)))
g<- function(x,a,b){
     out<-c()
     for (i in 1:length(x)){
       out[i]<-1- (integrate(model,0,1,x[i],a,b)$value / beta(a,b))
     }
     out
     }
x<- 10^seq(0,10,by=0.01)
y<- g(x,a=0.8,b=0.5)

This gives the error

Error in integrate(model, 0, 1, x[i], a, b) :   the integral is
probably divergent


Changing the relative or absolute tolerance solves this issue, but a certain
tolerance only works with a certain set of 'a' and 'b'.
For example, and abs.tol=10^-9 will make it work with a=0.8 and b=0.5 but
fail with a=0.3 and b=0.9.
We need this code to work for any "reasonable" value of 'a' and 'b' - as
seen by the shape of the distribution Beta(a,b).

We have tried using a different number of subdivisions without any luck.
The same integration in MATLAB works without any problem (using quad).

Anyone has an idea of why these problems occur and how to avoid them?

Many thanks.

	[[alternative HTML version deleted]]

