Hello,
We are trying to use R to simulate a model based on some 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]]
______________________________________________
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.
