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@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
______________________________________________
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.
