Dear R-experts,
I have 28 data points that I would like to fit with a non linear
broken-stick -- with three fitted parameters.
When I view trace -- and use the final values as lines on the graph of
data -- it looks pretty good.
Q1. Why am I getting singular convergeance?
Q2. Can you suggest another approach that may prove more satisfying?
I have read previous examples on nls and this sort of error on the list,
but can't find anything completely comparable.
Unfortunately the brocken stick function is a bit lengthy, but I have
cut down my example as much as possible.
Input
t <- c( 0.000000, 1.050000, 2.133333,
3.200000,4.266667,5.350000,6.416667,
7.483333,14.016667,15.133333,16.183333,17.300000,18.400000,19.466667
,20.550000,21.616667,22.716667,23.800000,30.350000,31.433333,32.533333
,33.616667,34.700000,35.783333,36.900000,38.000000,39.083333,40.166667)
seg_an <-c(
357.0466,360.5417,364.7510,368.8358,373.5498,377.2899,380.8858,385.9601
,420.2834,438.3807,453.9618,473.7764,493.0898,513.0759,531.1967,549.5310
,564.8920,584.8651,670.3014,674.2099,677.9492,680.5667,684.3941,688.2404
,690.7223,693.3406,697.9022,700.6606)
trans <- 13.38
trans2 <- 28.53
estCd <- 1975
estConst1 <- 0.00115689
estExch <- 0.00171680
Cb <- 330.1
Cmin <- 357
Ci <- 101000
R <- 0.001104768
A <- 16
V <- 8
rismod <- nls(seg_an ~ crv5(t, R, exch, trans, trans2, Cd, const1,
Cmin, Cb, A, V, Ci, Cmin.new), start = list(Cd = estCd, const1 =
estConst1,
exch = estExch), lower = c(100, 1e-05, 0),
upper = c(50000, 1, 0.0359), algorithm = "port", trace=TRUE)
plot(t, seg_an)
lines(t, crv5((t, R, 0.00182669, trans, trans2, 1976.6, 0.00115723,
Cmin, Cb, A, V, Ci, Cmin.new))
Output:
> rismod <- nls(seg_an ~ crv5(t, R, exch, trans, trans2, Cd, const1,
+ Cmin, Cb, A, V, Ci, Cmin.new), start = list(Cd = estCd, const1 =
estConst1,
+ exch = estExch), lower = c(100, 1e-05, 0),
+ upper = c(50000, 1, 0.0359), algorithm = "port", trace=TRUE)
0: 15.566264: 1975.00 0.00115689 0.00171680
1: 15.513993: 1983.53 0.00115234 0.00190479
2: 15.513949: 1979.94 0.00115487 0.00186440
3: 15.513949: 1976.35 0.00115739 0.00182390
4: 15.513946: 1978.15 0.00115613 0.00184450
5: 15.513946: 1976.35 0.00115739 0.00182418
6: 15.513946: 1977.25 0.00115676 0.00183432
7: 15.513946: 1976.80 0.00115707 0.00182908
8: 15.513946: 1976.58 0.00115723 0.00182669
9: 15.513946: 1976.58 0.00115723 0.00182669
10: 15.513946: 1976.58 0.00115723 0.00182669
Error in nls(seg_an ~ crv5(t, R, exch, trans, trans2, Cd, const1, Cmin,
:
Convergence failure: singular convergence (7)
crv5 function:
crv5 <- function(t, R, exch, trans, trans2, Cd, const1, Cmin=340,
Cb=325, A=16, V=8, Ci=93000, Cmin.new=NA) {
#function for broken stick nls regression, which joins to the three
#curves as a single function -- native emission rise and native+tracer
gas emission rise
pt1 <-
exp(-(t*(0+const1*A+exch))/V)*(((Ci*0+const1*Cd*A+Cb*exch)
*exp((t*(0+const1*A+exch))/V))/
(0+const1*A+exch)+((Cmin-Ci)*0+(const1*Cmin-const1*Cd)*A+(Cmin-Cb)*exch)
/(0+const1*A+exch))
if(is.na(Cmin.new)){
Cmin.new <-
exp(-(trans*(0+const1*A+exch))/V)*(((Ci*0+const1*Cd*A+Cb*exch)
*exp((trans*(0+const1*A+exch))/V))/
(0+const1*A+exch)+((Cmin-Ci)*0+(const1*Cmin-const1*Cd)*A+(Cmin-Cb)*exch)
/(0+const1*A+exch))
}
t.new <- t - trans
pt2 <-
exp(-(t.new*(R+const1*A+exch))/V)*(((Ci*R+const1*Cd*A+Cb*exch)
*exp((t.new*(R+const1*A+exch))/V))/
(R+const1*A+exch)+((Cmin.new-Ci)*R+(const1*Cmin.new-const1*Cd)*A+(Cmin.n
ew-Cb)*exch)/(R+const1*A+exch))
Cmin.new2 <-
exp(-((trans2-trans)*(R+const1*A+exch))/V)*(((Ci*R+const1*Cd*A+Cb*exch)
*exp(((trans2-trans)*(R+const1*A+exch))/V))/
(R+const1*A+exch)+((Cmin.new-Ci)*R+(const1*Cmin.new-const1*Cd)*A+(Cmin.n
ew-Cb)*exch)/(R+const1*A+exch))
t.new <- t - trans2
pt3 <-
exp(-(t.new*(0+const1*A+exch))/V)*(((Ci*0+const1*Cd*A+Cb*exch)
*exp((t.new*(0+const1*A+exch))/V))/
(0+const1*A+exch)+((Cmin.new2-Ci)*0+(const1*Cmin.new2-const1*Cd)*A+(Cmin
.new2-Cb)*exch)/(0+const1*A+exch))
res <- pt1*(t <=trans) + pt2 * ((t > trans) & (t <= trans2)) +
pt3 * (t > trans2)
return(res)
}
Thanks for your help.
Kind regards,
Matt Redding
We're behind the Bid!
GOLD COAST 2018 - XXI COMMONWEALTH GAMES CANDIDATE CITY
www.goldcoast2018bid.com
********************************DISCLAIMER**************...{{dropped:15}}
______________________________________________
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.
