Greetings List,
 
2 issues; likfit 
 
I'm having trouble with a simulation with 1,000 replications.  I'm trying to compare a series of models on differing simulated 2D surfaces.  The model below will have the following error at some iteration (so far the problems arises somewhere the 5th to 82th iteration).
 
The first section below is a portion of the code and the second section is the output and traceback().  I'm using variog in geoR to estimate sigmasq and phi to plug into the likfit function.  
 
Since the only values to change between iterations is the simulated 'y' and associated estimated sigmasq and phi, I assumed that the problem was in the sigmasq or phi.  Therefore, I manually changed the values on this iteration of the simulation to levels that seemed closer to an average (or use the ini.vals used in the variofit function).  When the values were set to these levels, likfit ran properly.  What I've found is that these estimated levels of sigmasq and phi that cause the problem are somewhat higher than the values that ran properly.  
 
To prevent this, I used the limits=pars.limits() in the call for variofit function.  I wanted to limit phi to be between 1.6 and 2.5 and limit sigmasq to be between 20 and 25.  When I set up pars.limits() that way I thought it should be, the on-screen output was correct but the variofit function read the upper limit for phi to be the upper limits for both phi AND sigmasq.
 
I would like to continue with my automated simuation.  Any suggestions are greatly appreciated,

Terry
 
 
 
 
The portion of the code is below (the y's were already simulated)
############
vario100 <- variog(coords=data.geo$coords, data=y, messages=FALSE)#estimates semivariogram
ini.vals <- c(21.1,2)
#limits=pars.limits(phi=c(2.5,25), sigmasq=c(20, 22.5)) ##I don't know why this doesn't work as expected so changed lower/upper bounds to make work
cressie <- variofit(vario100, ini=ini.vals, fix.nug=TRUE, wei="cressie",limits=pars.limits(phi=c(2.5,25), sigmasq=c(20, 22.5)), messages=FALSE)#estimates partial sill and range of semivariogram
#summary(cressie)#only as a test
#plot(vario100)#only as a test
#lines(cressie, lty=2)#only as a test
 
collVARIO<-cressie[2]$cov.pars#collects partial sill and range of semivariogram
geo.reml<-likfit(data.geo,coords=data.geo$coords,data=y,cov.model="gauss",trend=~s2+s3+s4+t+ts2+ts3+ts4,method.lik="REML",ini=c(cressie[2]$cov.pars[1],cressie[2]$cov.pars[2]), messages=FALSE)#geostatistical estimator using estimated partial sill and range
#geo.reml<-likfit(data.geo,coords=data.geo$coords,data=y,cov.model="gauss",trend=~s2+s3+s4+t+ts2+ts3+ts4,method.lik="REML",ini=c(21.94372,2.444176), messages=FALSE)#geostatistical estimator using estimated partial sill and range
###########

 
#################
#output from R
#######.
 
> geo.reml<-likfit(data.geo,coords=data.geo$coords,data=y,cov.model="gauss",trend=~s2+s3+s4+t+ts2+ts3+ts4,method.lik="REML",ini=c(21.94372,2.444176), messages=FALSE)#geostatistical estimator using estimated partial sill and range
trying another decomposition (svd)
Error in solve.default(v$varcov, xmat) : system is computationally singular: reciprocal condition number = 3.62888e-021
> traceback()
8: .Call("La_dgesv", a, b, tol, PACKAGE = "base")
7: solve.default(v$varcov, xmat)
6: solve(v$varcov, xmat)
5: solve(v$varcov, xmat)
4: fn(par, ...)
3: function (par) 
   fn(par, ...)(c(24.0717362048702, -1.56558793706907e-16))
2: optim(par = ini, fn = .negloglik.GRF, method = "L-BFGS-B", lower = lower.optim, 
       upper = upper.optim, fp = fixed.values, ip = ip, temp.list = temp.list, 
       ...)
1: likfit(data.geo, coords = data.geo$coords, data = y, cov.model = "gauss", 
       trend = ~s2 + s3 + s4 + t + ts2 + ts3 + ts4, method.lik = "REML", 
       ini = c(21.94372, 2.444176), messages = FALSE)
> 




Terry Griffin
Graduate Research Assistant
Agricultural Economics
Purdue University
765.586.2385
http://web.ics.purdue.edu/~twgriffi/

	[[alternative HTML version deleted]]


