[R] help with function

[Mike Nielsen]

Take the case of i==1.

Ct[i]<-1/Bq*Bt[i]*Cerr                                  # Assign Ct[1]
using Bt[1]
  Rt[i]<-Bt[i]/(a+b*Bt[i])                               # Assign
Rt[1] using Bt[1]
  Bt[i+2]<-(1-m)*Bt[i+1]+Rt[i]   *Rerr-Ct[i+1]  # Assign Bt[3] using
Bt[2] and Rt[1] and **Ct[2]**


You're reading Ct[i+1] before you ever assign it, hence NA.

OSISTM

Hope this helps,

Regards,

Mike



On 9/20/06, Guenther, Cameron <Cameron.Guenther at myfwc.com> wrote:
> Hello everyone,
>
> I have a function here that I wrote but doesn't seem to work quite
> right.  Attached is the code.  In the calib funcion under the for loop
> Bt[i+2]<-(1-m)*Bt[i+1]+Rt[i]*Rerr-Ct[i+1] returns NA's for everything
> after years 1983 and 1984.  However the code works when it reads
> Bt[i+2]<-(1-m)*Bt[i+1]+Rt[i]*Rerr-Ct[i].  I don't quite understand why
> since it should be calculating all of the necessary inputs prior to
> calculating Bt[i+2].  Any help would be greatly appreciated.
>
> Thanks
>
> #Model parameters
> B0<-75000000
> m<-0.3
> R0<-B0*m
> z<-0.8
> a<-B0/R0*(1-(z-0.2)/(0.8*z))
> b<-(z-0.2)/(0.8*z*R0)
> dat<-data.frame(years=seq(1983,2004),cobs=c(19032324,19032324,17531618,2
> 0533029,20298099,20793744,23519369,23131780,19922247,17274513,17034419,1
> 2448318,4551585,4226451,7183688,7407924,7538366,7336039,8869193,7902341,
> 6369089,6211886))
> stdr<-runif(100,0,0.5)
> stdc<-runif(100,0,0.5)
> BC<-runif(1000,0,100)
>
>
> #model calibration
>
> calib<-function(x){
>  v<-sample(stdr,1)
>  cr<-sample(stdc,1)
>  N<-rnorm(1)
>  Bq<-sample(BC,1)
>  Rerr<-exp(N*v-(v^2/2))
>  Cerr<-exp(N*cr-(cr^2/2))
>  Bt<-vector();Bt[1]=B0;Bt[2]=B0
>  Rt<-vector()
>  Ct<-vector()
>  for (i in 1:length(x$years)){
>   Ct[i]<-1/Bq*Bt[i]*Cerr
>   Rt[i]<-Bt[i]/(a+b*Bt[i])
>   Bt[i+2]<-(1-m)*Bt[i+1]+Rt[i]*Rerr-Ct[i+1]
>  }
>   out<-new.env()
>   out$yr<-x$years[1:length(x$years)]
>   out$Bt<-Bt[1:length(x$years)]
>   out$Rt<-Rt[1:length(x$years)]
>   out$Ct<-Ct[1:length(x$years)]
>   out$stdr<-v
>   out$stdc<-cr
>   out$Bq<-Bq
>   out$Rerr<-Rerr
>   out$Cerr<-Cerr
>   return(as.list(out))
>  }
>  test<-calib(dat)
>
>
> Cameron Guenther, Ph.D.
> Associate Research Scientist
> FWC/FWRI, Marine Fisheries Research
> 100 8th Avenue S.E.
> St. Petersburg, FL 33701
> (727)896-8626 Ext. 4305
> cameron.guenther at myfwc.com
>
> ______________________________________________
> 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
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Regards,

Mike Nielsen