[R] numerical integration problem
Sundar Dorai-Raj
sundar.dorai-raj at pdf.com
Thu Jun 29 13:38:55 CEST 2006
przeszczepan wrote:
> Hi,
>
> I have got problems integrating the following function using "integrate":
>
> lambdat<-function(t){
> tempT<-T[k,][!is.na(T[k,])]#available values from k-th row of matrix T
> tempJ<-J[k,][!is.na(J[k,])]
>
> hg<-length(tempT[tempT<=t & tempJ==0])#counts observations satisfing the conditions
> ag<-length(tempT[tempT<=t & tempJ==1])
>
> lambdaXY[hg+1,ag+1]#takes values from a 10x10 matrix
> }
>
> I keep receiving this message:
>
> 1: longer object length
> is not a multiple of shorter object length in: tempT <= t
> 2: longer object length
> is not a multiple of shorter object length in: tempT <= t & tempJ == 0
>
> What I suspect is that the "integrate" function submits the whole vector of points at which the integral is to be evaluated at once. For my function to be integrated it would rather have to be evaluated at each point after another in a loop of some kind.
>
You suspect correctly. Best to read ?integrate too.
> Can you think of a way to solve this problem without me having to write the integrating procedure from scratch (I have no idea about FORTRAN and this is what the "integrate" description refers to)?
>
Just put a "for"-loop in your function to iterate over t.
n <- length(t)
hg <- ag <- vector("numeric", n)
for(i in seq(n)) {
hg[i] <- length(tempT[tempT <= t[i] & tempJ == 0])
ag[i] <- length(tempT[tempT <= t[i] & tempJ == 1])
}
I doubt this will work because integrate is expecting a vector of
n=length(t) from lambdat. The last line of the function returns a nxn
matrix. Please submit data to run the function plus your call to
integrate in any subsequent postings.
HTH,
--sundar
> Thank you.
>
> Kind Regards,
> Lukasz Szczepanski
> Student
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
More information about the R-help
mailing list