[R] Integration problem: error in invoking an outside function

[Ravi Varadhan]

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.