[R] Newbie q. need some help understanding this code.
Kevin Wang
Kevin.Wang at maths.anu.edu.au
Thu Sep 16 02:34:49 CEST 2004
Hi,
sean kim wrote:
> thanks for any insights or other comments.
I would suggest that you run the code line by line, to see what it does
yourself. It is the best way to learn!
> library(stepfun)
This just loads the package "stepfun".
> lv <- function(N=1000,cvec=c(1,0.005,0.6),x=c(50,100))
> {
Declaration of a function named lv, where N, cvec and x are the
parameters/arguments with default values.
> m<-length(cvec)
m is a new variable, with in this case is just a number, the length of
the vector cvec. In this case, 3.
> n<-length(x)
Ditto.
> xmat<-matrix(nrow=N+1,ncol=n)
Generate a matrix with N + 1 (1001) rows and n columns.
> tvec<-vector("numeric",N)
> h<-vector("numeric",m)
Generating two vectors.
> t<-0
> xmat[1,]<-x
Assign x to the first row of the xmat matrix.
> for (i in 1:N) {
> h[1]<-cvec[1]*x[1]
The first element in cvec and x, multiply them together and the result
becomes the first element of h.
> h[2]<-cvec[2]*x[1]*x[2]
> h[3]<-cvec[3]*x[2]
Ditto.
> h0=sum(h)
Get the sum of h.
> tp<-rexp(1,h0)
Generating an exponential random number with a rate of h0.
> t<-t+tp
> u<-runif(1,0,1)
Generating a uniform random number, [0, 1]
> if ( u < h[1]/h0 ) {
> x[1] <- x[1]+1
If the uniform random number, u, is smaller than h[1]/h0 then do...
otherwise do the following...
> } else if ( u < (h[1]+h[2])/h0 ) {
> x[1] <- x[1]-1
> x[2] <- x[2]+1
> } else {
> x[2] <- x[2]-1
> }
> xmat[i+1,]<-x
Increment to the next row.
> tvec[i]<-t
> }
>
> list(stepfun(tvec,xmat[,1]),stepfun(tvec,xmat[,2]))
> }
Put things together into a list.
Again, check these line by line and you'll have a better understanding!
Kev
--
Ko-Kang Kevin Wang
PhD Student
Centre for Mathematics and its Applications
Building 27, Room 1004
Mathematical Sciences Institute (MSI)
Australian National University
Canberra, ACT 0200
Australia
Homepage: http://wwwmaths.anu.edu.au/~wangk/
Ph (W): +61-2-6125-2431
Ph (H): +61-2-6125-7407
Ph (M): +61-40-451-8301
More information about the R-help
mailing list