On 08-03-2012, at 14:50, Journals wrote: > I wish to fit a dynamical model in R and I am running in a problem that > requires some of your wisdom to solve. For SAS users I am searching for > the equivalent of the */retain/ *statement. > > For people that want to read complicated explanations to help me: > > I have a system of two equations written as difference equations here. > > To boil it down. I have a dataframe with three variables y, X1, X2 which > are measured. Now, I want to estimate a model which says y= > yhat+epsillon as in standard nonlinear regression. Where yhat is the > predicted value of y. > > yhat can be calculated as follows: > > # This code illustrates how I could simulate the expected values of y if > I knew the values of the parameters tau and b. > # but in reality I would like to estimate them. > # code is for illustration of the principles and is not meant to be > functional!! > yaux[1]<-0 > b<- a_number # b would have to become estimated by nls or nlme > tau<- another_number # tau would also be estimated in nls or nlme > for (t in 2:1000) > { > yaux[t+1] <- yaux[t] + (X1-yaux[t])/tau > yhat[t+1] <- yaux[t+1]*X2[t+1]/(X2[t+1]+b) > } >
This can't be correct. When t==2 you are referencing yaux[2] which has no value. Secondly, you haven't indexed X1 in the equation for yaux[t+1]. Is it X1[t]? Besides that you don't have to calculate yaux that way. You can use filter. As in # Initial value yaux is 0 yaux <- as.vector(filter(a*X1,c(1-a),method="recursive")) where a = 1/tau and assuming you meant X1[t]. And then you don't need to calculate yhat recursively. You can do it with yhat <- yaux * X2/(X2+b) > Now, my problem is that I do not know the values of /tau/ and /b/ and I > would like to estimated them by non-linear regression. This is easy in > the case of /b/ but for /tau /nls (or nlme etc) would have to remember > the value of /yaux /for the previous observation and I did not find any > syntactical mean to do that. > You could do a non linear regression with y.rhs <- function(a,b,X1,X2) as.vector(filter(a*X1,c(1-a),method="recursive")) * X2/(X2+b) nls(y ~ y.rhs(a,b,X1,X2), start=list(a=a,b=b)) where tau=1/a Here is an example with random data N <- 500 a <- .4 b <- 1 # generate some random data set.seed(1) X1 <- 5*round(runif(N),2) X2 <- round(runif(N),2) yaux <- as.vector(filter(a*X1,c(1-a),method="recursive")) yhat <- yaux * X2/(X2+b) eps <- runif(N)/2 eps <- eps - mean(eps) y <- yhat + eps sum((y-yhat)^2) # starting values a <- .6 b <- 2 y.rhs <- function(a,b,X1,X2) as.vector(filter(a*X1,c(1-a),method="recursive")) * X2/(X2+b) nls(y ~ y.rhs(a,b,X1,X2), start=list(a=a,b=b)) Berend ______________________________________________ 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.