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

[David Winsemius]

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-matrix

Not 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