Frankvb
2010-Jun-15 13:05 UTC
[R] Integration problem: error in invoking an outside function
Dear all, Currently I am trying to integrate a function which depends on four variables, two of which are given, one is given in the integrate function, so there is one variable to integrate on. The code is as follows: Pmatrix = function(th) { P = matrix(nrow=6, ncol=6, data=0) P[1,1] = P[2,1]=P[3,2]=P[4,3]=P[5,4]=P[6,5]= exp(-th) P[,6] = 1-exp(-th) return(P)} lim.verd = function(matrix) { et = matrix(nrow=1, ncol=dim(matrix)[2], data=1) E = matrix(nrow=dim(matrix)[1], ncol=dim(matrix)[2], data=1) mat = diag(dim(matrix)[1]) - matrix + E inverse.mat = solve(mat) pi = et %*% inverse.mat return(pi)} a.hat = 0.8888 lambda.hat = 0.1474 int = function(theta, s, a, lambda) { a = a.hat lambda = lambda.hat f.dist = gamma(a)^(-1) * a^a * theta^(a-1) * exp(-a*theta) # numeric pi = lim.verd(PM((lambda*theta))) # a 1x6-matrix return(pi[1,s+1]*f.dist)} integrate(int,lower=0.0001,upper=10,s=2) If I try int(0.1,2) a get a numerical result. However, once I use integrate, I get the following.> integrate(int,lower=0.0001,upper=10,s=2)Error in P[1, 1] = exp(-th) : number of items to replace is not a multiple of replacement length I've searched the internet, but I haven't been able to find a solution to this yet. Does anyone have a clue? Thanks in advance! Frank van Berkum -- View this message in context: http://r.789695.n4.nabble.com/Integration-problem-error-in-invoking-an-outside-function-tp2255862p2255862.html Sent from the R help mailing list archive at Nabble.com.
Ravi Varadhan
2010-Jun-15 13:49 UTC
[R] Integration problem: error in invoking an outside function
Frank, The trouble is that your integrand function `int' is not vectorized. You can use `Vectorize' or `sapply' to rectify this. int.vec = function(theta, s, a, lambda) { sapply(theta, function(t,s,a,lambda) int(t,s,a,lambda), s=s, a=a, lambda=lambda) } integrate(int.vec,lower=0.0001,upper=10,s=2) When I do this I get the following answer:> integrate(int.vec,lower=0.0001,upper=10,s=2)0.06224828 with absolute error < 9.4e-05>Note: The code that you posted is incorrect and not reproducible. pi = lim.verd(PM((lambda*theta))) # a 1x6-matrix The function `PM' is undefined. This should have been: pi = lim.verd(Pmatrix((lambda*theta))) # a 1x6-matrix Hope this helps, Ravi. -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Frankvb Sent: Tuesday, June 15, 2010 9:06 AM To: r-help at r-project.org Subject: [R] Integration problem: error in invoking an outside function Dear all, Currently I am trying to integrate a function which depends on four variables, two of which are given, one is given in the integrate function, so there is one variable to integrate on. The code is as follows: Pmatrix = function(th) { P = matrix(nrow=6, ncol=6, data=0) P[1,1] = P[2,1]=P[3,2]=P[4,3]=P[5,4]=P[6,5]= exp(-th) P[,6] = 1-exp(-th) return(P)} lim.verd = function(matrix) { et = matrix(nrow=1, ncol=dim(matrix)[2], data=1) E = matrix(nrow=dim(matrix)[1], ncol=dim(matrix)[2], data=1) mat = diag(dim(matrix)[1]) - matrix + E inverse.mat = solve(mat) pi = et %*% inverse.mat return(pi)} a.hat = 0.8888 lambda.hat = 0.1474 int = function(theta, s, a, lambda) { a = a.hat lambda = lambda.hat f.dist = gamma(a)^(-1) * a^a * theta^(a-1) * exp(-a*theta) # numeric pi = lim.verd(PM((lambda*theta))) # a 1x6-matrix return(pi[1,s+1]*f.dist)} integrate(int,lower=0.0001,upper=10,s=2) If I try int(0.1,2) a get a numerical result. However, once I use integrate, I get the following.> integrate(int,lower=0.0001,upper=10,s=2)Error in P[1, 1] = exp(-th) : number of items to replace is not a multiple of replacement length I've searched the internet, but I haven't been able to find a solution to this yet. Does anyone have a clue? Thanks in advance! Frank van Berkum -- View this message in context: http://r.789695.n4.nabble.com/Integration-problem-error-in-invoking-an-outsi de-function-tp2255862p2255862.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ R-help at 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.
David Winsemius
2010-Jun-15 14:17 UTC
[R] Integration problem: error in invoking an outside function
On Jun 15, 2010, at 9:05 AM, Frankvb wrote:> > Dear all, > > Currently I am trying to integrate a function which depends on four > variables, two of which are given, one is given in the integrate > function, > so there is one variable to integrate on. > > The code is as follows: > Pmatrix > function(th) { > P = matrix(nrow=6, ncol=6, data=0) > P[1,1] = P[2,1]=P[3,2]=P[4,3]=P[5,4]=P[6,5]= exp(-th) > P[,6] = 1-exp(-th) > return(P)} > > lim.verd > function(matrix) { > et = matrix(nrow=1, ncol=dim(matrix)[2], data=1) > E = matrix(nrow=dim(matrix)[1], ncol=dim(matrix)[2], data=1) > mat = diag(dim(matrix)[1]) - matrix + E > inverse.mat = solve(mat) > pi = et %*% inverse.mat > return(pi)} > > a.hat = 0.8888 > lambda.hat = 0.1474 > int > function(theta, s, a, lambda) { > a = a.hat > lambda = lambda.hat > f.dist = gamma(a)^(-1) * a^a * theta^(a-1) * exp(-a*theta) # numeric > pi = lim.verd(PM((lambda*theta))) # a 1x6-matrixNot a good idea to redefine pi. It's been tried before with laughable results.> return(pi[1,s+1]*f.dist)} > integrate(int,lower=0.0001,upper=10,s=2) > If I try int(0.1,2) a get a numerical result. However, once I use > integrate, > I get the following. >> integrate(int,lower=0.0001,upper=10,s=2) > Error in P[1, 1] = exp(-th) : > number of items to replace is not a multiple of replacement length > I've searched the internet, but I haven't been able to find a > solution to > this yet. Does anyone have a clue?The integrate function generates a sequence, ie. a vector, passed to the first argument, theta in this case, which is then multiplied by a scalar, lambda, before passing it as an argument to PM (which I assume is what you posted as Pmatrix) where it is given to exp() and then an attempt is made to assign to P[6,5] which fails, throwing the error. Learn some debugging strategies. -- David.> > Thanks in advance! > > Frank van Berkum > > -- > View this message in context: http://r.789695.n4.nabble.com/Integration-problem-error-in-invoking-an-outside-function-tp2255862p2255862.html > Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > R-help at 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.David Winsemius, MD West Hartford, CT