R help - Feb 2014 - Parameter fitting by using minpack.lm package

Hi!

I am trying to fit 2-dimensional data (global data with 2-degree
resolution) using nls.lm function of minpack.lm package. I am successful to
fit the transient data in a single grid point but not able to fit 2-D data.
The error message is listed below. I am just wondering if anybody has used
this routine to fit 2-D/3-D data?

Error in `colnames<-`(`*tmp*`, value = c("a", "b",
"c")) :
  length of 'dimnames' [2] not equal to array extent

It seems like the error is generated after the iteration. I googled the
error but did not find anything significant.

It.    0, RSS =        nan, Par. =          0          0          0
It.    1, RSS =       37.5, Par. =          0          0          0

Here is the code I used:

require(graphics)
library(ncdf)
library(minpack.lm)

data <-
"E:/R/Optimize_FIN.Params/output.with.zwt_actual.base.jul1993.nc"
dataModel <- open.ncdf(data)
data.obs <-
"E:/R/Optimize_FIN.Params/Prigent_fraction_jul1993_conserved.nc"
dataObs <- open.ncdf(data.obs)
data.params <- "E:/R/Optimize_FIN.Params/
surfdata_1.9x2.5_simyr1850_c120622.nc"
dataParams <- open.ncdf(data.params)

# get latitudes and longitudes
lat <- get.var.ncdf(dataModel,"lat",verbose=F)
nlat <- dim(lat)
lon <- get.var.ncdf(dataModel,"lon")
nlon <- dim(lon)
dimlat <- dim.def.ncdf("lat","degrees_north",lat)
dimlon <- dim.def.ncdf("lon","degrees_east",lon)
print(dimlat)

# get time dimension for model data
t <- get.var.ncdf(dataModel,"time")
nt <- dim(t)
tunits <- att.get.ncdf(dataModel,"time","units")

#get the variables from the netcdf file
zwt0 <- get.var.ncdf(dataParams,"ZWT0")
f0 <- get.var.ncdf(dataParams,"F0")
p3 <- get.var.ncdf(dataParams,"P3")
zwt <- get.var.ncdf(dataModel,"ZWT_ACTUAL")
qr <- get.var.ncdf(dataModel,"QOVER_LAG")
obs <- get.var.ncdf(dataObs,"frac")
close.ncdf(dataModel)
close.ncdf(dataParams)
close.ncdf(dataObs)

#change NA to zero in data
qr[is.na(qr)] <- 0
zwt[is.na(zwt)]<- 0

for (i in 1:144)
{
  for (j in 1:96)
  {
    f.fin <- function(params,data) {
      fin.mod <- expression(params$a*exp(-(data$a/params$b)) +
params$c*data$b)
      eval(fin.mod)
    }
    p.fin <- list(a=f0, b=zwt0, c=p3)
    d.fin <- list(a=zwt, b=qr)
    model <- f.fin(p.fin, d.fin)

    residFun <- function(p, observed, data) {
      r.residual <-expression(observed - f.fin(p,data))
      eval(r.residual)
    }
    parStart <- list(a=f0,b=zwt0,c=p3)
    #perform the optimization
    nls.out <- nls.lm(par=parStart, fn = residFun, observed = obs,
                      data = d.fin, control = nls.lm.control(maxiter
100,nprint=1))
  }
}
summary(nls.out)

Any help would be greatly appreciated. Thank you.

Best,
Rajendra Paudel
Cornell University

	[[alternative HTML version deleted]]