Delphine Pessoa
2011-Aug-03 15:07 UTC
[R] Case-by-case tolerance needed for successful integrate()
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]]
Delphine Pessoa
2011-Aug-17 18:55 UTC
[R] Case-by-case tolerance needed for successful integrate()
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]]