Hi. I am attempting to run my code to get estimates for a nonlinear
model. Unfortunately, when I run the code, I keep getting the following
errors (they switch back and forth depending on when I run it):
Error in nlsModel(formula, mf, start, wts) :
singular gradient matrix at initial parameter estimates
Error in nls(ywt_vib ~ weightfunc_vib(time_vib, Wmax_star, tau_star,
gamma_star, :
step factor 0.000488281 reduced below 'minFactor' of 0.000976563
I'll put my code at the bottom of this email. Any help you can offer
regarding this error would be greatly appreciated. Thank you so much!!
Trena Phipps
-------------------------------------------------------------------------------------------
vib <- matrix(scan(),ncol=3,byrow=TRUE)
1 30.0 41.1
1 45.0 51.5
1 60.0 54.0
1 75.0 56.2
1 90.0 59.1
1 120.0 59.3
time_vib <- vib[,2]
diss_vib <- vib[,3]
outfile <- "work"
cat("WORK","\n",file=outfile,"\n",append=TRUE)
n_vib <- length(time_vib)
p_vib <- 3
cmax <- 30
WLSiter <- 500
beta0_vib <- list(Wmax_star= 30, tau_star=15, gamma_star=.8)
t0 <- 20
tol <- 10^-8
meanfunc_vib <- function(time_vib,Wmax_star,tau_star,gamma_star){
diff<-time_vib-t0
f_vib<-Wmax_star*(1-exp(-log(2)*((diff/tau_star)^gamma_star)))
f_vib
# meanfunc_vib_deriv <- function(time_vib,Wmax_star,tau_star,gamma_star)
meangrad_vib <- array(0,c(n_vib,p_vib))
meangrad_vib[,1] <- 1-exp(-log(2)*(diff/tau_star)^gamma_star)
meangrad_vib[,2] <-
-(Wmax_star*log(2)*((diff/tau_star)^gamma_star)*exp(-log(2)*((diff/tau_star)^gamma_star)))/tau_star
meangrad_vib[,3] <-
Wmax_star*log(2)*((diff/tau_star)^gamma_star)*log(diff/tau_star)*exp(-log(2)*((diff/tau_star)^gamma_star))
attr(f_vib,"gradient") <- meangrad_vib
f_vib
}
unweightfunc_vib <- function(time_vib,Wmax_star,tau_star,gamma_star){
diff<-time_vib-t0
f_vib<-Wmax_star*(1-exp(-log(2)*(diff/tau_star)^gamma_star))
f_vib
}
weightfunc_vib <- function(time_vib,Wmax_star,tau_star,gamma_star,wt_vib){
diff<-time_vib-t0
f_vib<-Wmax_star*(1-exp(-log(2)*(diff/tau_star)^gamma_star))
w12_vib <- sqrt(wt_vib)
weightf_vib <- f_vib*w12_vib
weightgrad_vib <-
array(0,c(n_vib,p_vib),list(NULL,c("Wmax_star","tau_star", "gamma_star")))
weightgrad_vib[,"Wmax_star"] <-
w12_vib*(1-exp(-log(2)*(diff/tau_star)^gamma_star))
weightgrad_vib[,"tau_star"] <-
-w12_vib*((Wmax_star*log(2)*((diff/tau_star)^gamma_star)*exp(-log(2)*(diff/tau_star)^gamma_star))/tau_star)
weightgrad_vib[,"gamma_star"] <-
w12_vib*(Wmax_star*log(2)*((diff/tau_star)^gamma_star)*log(diff/tau_star)*exp(-log(2)*(diff/tau_star)^gamma_star))
attr(weightf_vib,"gradient") <- weightgrad_vib
weightf_vib
}
trkfunc_vib <- function(resid_vib,mudot_vib,mu_vib,theta_vib){
trk_vib <- resid_vib*((mudot_vib/mu_vib)**theta_vib)
trkgrad_vib <- array(0,c(length(mu_vib),1),list(,NULL,c("theta_vib")))
trkgrad_vib[,"theta_vib"] <- trk_vib*log(mudot_vib/mu_vib)
attr(trk_vib,"gradient") <- trkgrad_vib # analytic derivative
trk_vib
}
diss_vib_update<-diss_vib-30.50973
vib_dat<-data.frame(time_vib,diss_vib_update)
ols.fit_vib <- nls(diss_vib_update ~
meanfunc_vib(time_vib,Wmax_star,tau_star,gamma_star),vib_dat,beta0_vib)
bols_vib <- coef(ols.fit_vib)
mu_vib <- unweightfunc_vib(time_vib,bols_vib[1],bols_vib[2],bols_vib[3])
resid_vib <- diss_vib_update-mu_vib
sigols_vib <- sqrt(sum(resid_vib*resid_vib)/(n_vib-p_vib))
cat("FIT BY GLS",file=outfile,"\n",
append=F)
cat("OLS estimate - VIB ",round(bols_vib,6),file=outfile,"\n","\n",append=T)
cat("Estimate of assumed constant SD of y
",round(sigols_vib,6),file=outfile,
"\n","\n","\n",append=T)
bgls_vib <- bols_vib
diss_vib_update2<-diss_vib-30.57897
vib_dat<-data.frame(time_vib,diss_vib_update2)
for (k in 1:cmax){
mu_vib <- unweightfunc_vib(time_vib,bgls_vib[1],bgls_vib[2],bgls_vib[3])
mudot_vib <- prod(mu_vib)**(1/n_vib)
mudot_vib <-rep(mudot_vib,n_vib)
resid_vib <- diss_vib_update2-mu_vib
dummy_vib <- rep(0,length(time_vib))
pldat_vib <- data.frame(resid_vib,mu_vib,mudot_vib)
theta_pl_est_vib <- nls(dummy_vib ~
trkfunc_vib(resid_vib,mudot_vib,mu_vib,theta_vib),pldat_vib,list(theta_vib=.8),nls.control(maxiter=300))
theta_vib <- coef(theta_pl_est_vib)
mu_vib <- unweightfunc_vib(time_vib,bgls_vib[1],bgls_vib[2],bgls_vib[3])
wt_vib <- 1/(mu_vib^(2*theta_vib))
ywt_vib <- diss_vib_update2*sqrt(wt_vib)
thedat2_vib <- data.frame(time_vib,ywt_vib,wt_vib)
gls_est_vib <- nls(ywt_vib ~
weightfunc_vib(time_vib,Wmax_star,tau_star,gamma_star,wt_vib),thedat2_vib,beta0_vib,nls.control(maxiter=200))
bgls_vib <- coef(gls_est_vib)
cat("VIB - Iteration ",k,"\n",file=outfile,append=TRUE)
cat("Estimate of theta - VIB
",round(theta_vib,8),"\n",file=outfile,append=TRUE)
cat("GLS estimate of beta - VIB
",round(bgls_vib,8),"\n","\n",file=outfile,append=TRUE)
}
mu_vib <- unweightfunc_vib(time_vib,bgls_vib[1],bgls_vib[2],bgls_vib[3])
resid_vib <- diss_vib_update2-mu_vib
g_vib <- mu_vib^theta_vib
sigma2_vib <- sum((resid_vib/g_vib)**2)/(n_vib-p_vib)
sigma_vib <- sqrt(sigma2_vib)
cat("Final estimate of sigma^2 -
VIB",round(sigma2_vib,10),"\n","\n",file=outfile,append=TRUE)
sink(outfile,append=TRUE)
print(summary(gls_est_vib))
sink()
______________________________________________
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.
