Hello there,
I am trying to fit an exponential fit using Least squares to some data. #data x <- c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) y <- c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, -0.000006, -0.004626, -0.004626, -0.004626, -0.004626) sub <- data.frame(x,y) #If model is y = a*exp(-x) + b then fit <- nls(y ~ a*exp(-x) + b, data = sub, start = list(a = 0, b = 0), trace = TRUE) This works well. However, if I want to fit the model : y = a*exp(-mx)+c then I try - fit <- nls(y ~ a*exp(-m*x) + b, data = sub, start = list(a = 0, b = 0, m= 0), trace = TRUE) It fails and I get the following error - Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates Any suggestions how I can fix this? Also next I want to try to fit a sum of 2 exponentials to this data. So the new model would be y = a*exp[(-m1+ m2)*x]+c . Any suggestion how I can do this... Any help would be most appreciated. Thanks in advance. Diviya [[alternative HTML version deleted]] ______________________________________________ [hidden email] 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. |
On Jun 12, 2011, at 18:57 , Diviya Smith wrote: > Hello there, > > I am trying to fit an exponential fit using Least squares to some data. > > #data > x <- c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) > y <- c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, > -0.000006, -0.004626, -0.004626, -0.004626, -0.004626) > > sub <- data.frame(x,y) > > #If model is y = a*exp(-x) + b then > fit <- nls(y ~ a*exp(-x) + b, data = sub, start = list(a = 0, b = 0), trace > = TRUE) > > This works well. However, if I want to fit the model : y = a*exp(-mx)+c then > I try - > fit <- nls(y ~ a*exp(-m*x) + b, data = sub, start = list(a = 0, b = 0, m= > 0), trace = TRUE) > > It fails and I get the following error - > Error in nlsModel(formula, mf, start, wts) : > singular gradient matrix at initial parameter estimates If a==0, then a*exp(-m*x) does not depend on m. So don't use a=0 as initial value. > > Any suggestions how I can fix this? Also next I want to try to fit a sum of > 2 exponentials to this data. So the new model would be y = a*exp[(-m1+ > m2)*x]+c . That's not a sum of exponentials. Did you mean a*(exp(-m1*x) + exp(-m2*x)) + c? Anyways, same procedure with more parameters. Just beware the fundamental exchangeability of m1 and m2, so don't initialize them to the same value. > Any suggestion how I can do this... Any help would be most > appreciated. Thanks in advance. -- Peter Dalgaard Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: [hidden email] Priv: [hidden email] ______________________________________________ [hidden email] 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. |
On Sun, 12 Jun 2011, peter dalgaard wrote:
> > On Jun 12, 2011, at 18:57 , Diviya Smith wrote: > >> Hello there, >> >> I am trying to fit an exponential fit using Least squares to some data. >> >> #data >> x <- c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) >> y <- c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, >> -0.000006, -0.004626, -0.004626, -0.004626, -0.004626) >> >> sub <- data.frame(x,y) >> >> #If model is y = a*exp(-x) + b then >> fit <- nls(y ~ a*exp(-x) + b, data = sub, start = list(a = 0, b = 0), trace >> = TRUE) >> >> This works well. However, if I want to fit the model : y = a*exp(-mx)+c then >> I try - >> fit <- nls(y ~ a*exp(-m*x) + b, data = sub, start = list(a = 0, b = 0, m= >> 0), trace = TRUE) >> >> It fails and I get the following error - >> Error in nlsModel(formula, mf, start, wts) : >> singular gradient matrix at initial parameter estimates > > > If a==0, then a*exp(-m*x) does not depend on m. So don't use a=0 as initial value. > >> >> Any suggestions how I can fix this? Also next I want to try to fit a sum of >> 2 exponentials to this data. So the new model would be y = a*exp[(-m1+ >> m2)*x]+c . > > That's not a sum of exponentials. Did you mean a*(exp(-m1*x) + > exp(-m2*x)) + c? OK, that makes more sense. Also, scaling of the variables may help. Something like this could work: ## scaled data d <- data.frame(x = c(1, 1:10 * 10), y = 100 * c( 0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, -0.000006, -0.004626, -0.004626, -0.004626, -0.004626)) ## model fits n1 <- nls(y ~ a*exp(-m * x) + b, data = d, start = list(a = 1, b = 0, m = 0.1)) n2 <- nls(y ~ a * (exp(-m1 * x) + exp(-m2 * x)) + b, data = d, start = list(a = 1, b = 0, m1 = 0.1, m2 = 0.5)) ## visualization plot(y ~ x, data = d) lines(d$x, fitted(n1), col = 2) lines(d$x, fitted(n2), col = 4) ## ANOVA anova(n1, n2) ## LR test library("lmtest") lrtest(n1, n2) which both seem to indicate that n1 is sufficient. hth, Z > Anyways, same procedure with more parameters. Just > beware the fundamental exchangeability of m1 and m2, so don't initialize > them to the same value. > >> Any suggestion how I can do this... Any help would be most >> appreciated. Thanks in advance. > > -- > Peter Dalgaard > Center for Statistics, Copenhagen Business School > Solbjerg Plads 3, 2000 Frederiksberg, Denmark > Phone: (+45)38153501 > Email: [hidden email] Priv: [hidden email] > > ______________________________________________ > [hidden email] 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. > ______________________________________________ [hidden email] 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. |
In reply to this post by Diviya Smith
Hi:
If you use RStudio, then you can use its manipulate package to figure out starting values for the model visually through sliders. Walmes Zeviani posted the template of how to do this last week; I've just adapted it for the models under consideration. It's interesting to play with this because it shows how strongly some of the parameters are correlated with one another...and it's fun :) Thanks to Walmes for the excellent tutorial example. Firstly, x and y (or the inputs in general) need to be in the global environment as vectors of the same length. [RStudio is a GUI, so I'm assuming that you run things from its command line.] Load the manipulate package (which comes with RStudio) and use the manipulate() function to create a plot of the data and fit a curve to it. The sliders adjust the parameter values and you play with them until the curve fits closely to the data. The range of the sliders should match the range of plausible values you think they might take. The output of the manipulate() function is a vector of starting values that you can pass into nls(), which is done by closing the window containing the sliders. # Model 1: x <- c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) y <- c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, -0.000006, -0.004626, -0.004626, -0.004626, -0.004626) plot(x, y) # Initial plot of the data start <- list() # Initialize an empty list for the starting values library(manipulate) manipulate( { plot(x, y) k <- kk; b0 <- b00; b1 <- b10 curve(k*exp(-b1*x) + b0, add=TRUE) start <<- list(k=k, b0=b0, b1=b1) }, kk=slider(0.01, 0.2, step = 0.01, initial = 0.1), b10=slider(0.01, 0.4, step = 0.01, initial = 0.3), b00=slider(-0.01, -0.001, step=0.001,initial= -0.005)) # When done, start() is a list of named parameters in # the global environment # Model fit: fit1 <- nls(y ~ k*exp(-b1*x) + b0, start = start) summary(fit1) ### Model 2: [following Peter Dalgaard's suggested model] ### Use the estimates from fit1 and shrink b1 to anticipate ### potential effect of b2; the initial estimate in the slider for ### b1 will be too big, so it needs to be shrunk (a lot :) start <- list() manipulate( { plot(x, y) k <- kk; b0 <- b00; b1 <- b10; b2 <- b20 curve(k*(exp(-b1*x) + exp(-b2*x)) + b0, add=TRUE) start <<- list(k=k, b0=b0, b1=b1, b2 = b2) }, kk=slider(0.01, 0.1, step = 0.01, initial = 0.04), b10=slider(0.01, 0.4, step = 0.01, initial = 0.3), b20 = slider(0.01, 0.1, step = 0.005, initial = 0.05), b00=slider(-0.01, -0.001, step=0.001,initial= -0.004)) fit2 <- nls(y ~ k*(exp(-b1*x) + exp(-b2*x)) + b0, start = start) summary(fit2) anova(fit1, fit2) HTH, Dennis On Sun, Jun 12, 2011 at 9:57 AM, Diviya Smith <[hidden email]> wrote: > Hello there, > > I am trying to fit an exponential fit using Least squares to some data. > > #data > x <- c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) > y <- c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, > -0.000006, -0.004626, -0.004626, -0.004626, -0.004626) > > sub <- data.frame(x,y) > > #If model is y = a*exp(-x) + b then > fit <- nls(y ~ a*exp(-x) + b, data = sub, start = list(a = 0, b = 0), trace > = TRUE) > > This works well. However, if I want to fit the model : y = a*exp(-mx)+c then > I try - > fit <- nls(y ~ a*exp(-m*x) + b, data = sub, start = list(a = 0, b = 0, m= > 0), trace = TRUE) > > It fails and I get the following error - > Error in nlsModel(formula, mf, start, wts) : > singular gradient matrix at initial parameter estimates > > Any suggestions how I can fix this? Also next I want to try to fit a sum of > 2 exponentials to this data. So the new model would be y = a*exp[(-m1+ > m2)*x]+c . Any suggestion how I can do this... Any help would be most > appreciated. Thanks in advance. > > Diviya > > [[alternative HTML version deleted]] > > ______________________________________________ > [hidden email] 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. > ______________________________________________ [hidden email] 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. |
Thank you for very informative answers. This is great help. Peter, thanks
for correcting the model. That was exactly what I meant - apologies for the typo. I was able to run the analysis on some simulated data and the subset of data that I had posted earlier. However, when I apply the analysis to full dataset, I usually get the following errors- fm2 <- nls(y ~ b + a * (exp(m1 * -x) + exp(m2 * -x)), start = list(a = 1, b = 0, m1 = 0.1, m2 = 0.5), trace = TRUE) 84.42045 : 1.0 0.0 0.1 0.5 0.1593364 : 0.072399953 -0.000226868 0.101135423 0.571436093 0.09032837 : 7.897930e-02 7.905359e-06 1.271435e-01 1.872206e+00 0.06834495 : 0.0808842841 0.0005815865 0.1725302473 4.9048478439 Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model OR Error in nls(y ~ b + a * (exp(m1 * -x) + exp(m2 * -x)), start = list(a = 1, : number of iterations exceeded maximum of 50 Is there a way to optimize the start values? I was unable to try the RStudio solution as RStudio keeps crashing on my computer. Also is there any way to set the start value in NLS - I mean if I want to fit the model to only a subset of data starting at x = 10, how can I specify that? Thanks, Diviya On Sun, Jun 12, 2011 at 7:19 PM, Dennis Murphy <[hidden email]> wrote: > Hi: > > If you use RStudio, then you can use its manipulate package to figure > out starting values for the model visually through sliders. Walmes > Zeviani posted the template of how to do this last week; I've just > adapted it for the models under consideration. It's interesting to > play with this because it shows how strongly some of the parameters > are correlated with one another...and it's fun :) Thanks to Walmes for > the excellent tutorial example. > > Firstly, x and y (or the inputs in general) need to be in the global > environment as vectors of the same length. > [RStudio is a GUI, so I'm assuming that you run things from its > command line.] Load the manipulate package (which comes with RStudio) > and use the manipulate() function to create a plot of the data and fit > a curve to it. The sliders adjust the parameter values and you play > with them until the curve fits closely to the data. The range of the > sliders should match the range of plausible values you think they > might take. The output of the manipulate() function is a vector of > starting values that you can pass into nls(), which is done by closing > the window containing the sliders. > > # Model 1: > > x <- c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) > y <- c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, > -0.000006, -0.004626, -0.004626, -0.004626, -0.004626) > plot(x, y) # Initial plot of the data > start <- list() # Initialize an empty list for the starting values > > library(manipulate) > manipulate( > { > plot(x, y) > k <- kk; b0 <- b00; b1 <- b10 > curve(k*exp(-b1*x) + b0, add=TRUE) > start <<- list(k=k, b0=b0, b1=b1) > }, > kk=slider(0.01, 0.2, step = 0.01, initial = 0.1), > b10=slider(0.01, 0.4, step = 0.01, initial = 0.3), > b00=slider(-0.01, -0.001, step=0.001,initial= -0.005)) > > # When done, start() is a list of named parameters in > # the global environment > # Model fit: > fit1 <- nls(y ~ k*exp(-b1*x) + b0, start = start) > summary(fit1) > > ### Model 2: [following Peter Dalgaard's suggested model] > ### Use the estimates from fit1 and shrink b1 to anticipate > ### potential effect of b2; the initial estimate in the slider for > ### b1 will be too big, so it needs to be shrunk (a lot :) > start <- list() > manipulate( > { > plot(x, y) > k <- kk; b0 <- b00; b1 <- b10; b2 <- b20 > curve(k*(exp(-b1*x) + exp(-b2*x)) + b0, add=TRUE) > start <<- list(k=k, b0=b0, b1=b1, b2 = b2) > }, > kk=slider(0.01, 0.1, step = 0.01, initial = 0.04), > b10=slider(0.01, 0.4, step = 0.01, initial = 0.3), > b20 = slider(0.01, 0.1, step = 0.005, initial = 0.05), > b00=slider(-0.01, -0.001, step=0.001,initial= -0.004)) > > fit2 <- nls(y ~ k*(exp(-b1*x) + exp(-b2*x)) + b0, start = start) > summary(fit2) > > anova(fit1, fit2) > > > HTH, > Dennis > > On Sun, Jun 12, 2011 at 9:57 AM, Diviya Smith <[hidden email]> > wrote: > > Hello there, > > > > I am trying to fit an exponential fit using Least squares to some data. > > > > #data > > x <- c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) > > y <- c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, > > -0.000006, -0.004626, -0.004626, -0.004626, -0.004626) > > > > sub <- data.frame(x,y) > > > > #If model is y = a*exp(-x) + b then > > fit <- nls(y ~ a*exp(-x) + b, data = sub, start = list(a = 0, b = 0), > trace > > = TRUE) > > > > This works well. However, if I want to fit the model : y = a*exp(-mx)+c > then > > I try - > > fit <- nls(y ~ a*exp(-m*x) + b, data = sub, start = list(a = 0, b = 0, m= > > 0), trace = TRUE) > > > > It fails and I get the following error - > > Error in nlsModel(formula, mf, start, wts) : > > singular gradient matrix at initial parameter estimates > > > > Any suggestions how I can fix this? Also next I want to try to fit a sum > of > > 2 exponentials to this data. So the new model would be y = a*exp[(-m1+ > > m2)*x]+c . Any suggestion how I can do this... Any help would be most > > appreciated. Thanks in advance. > > > > Diviya > > > > [[alternative HTML version deleted]] > > > > ______________________________________________ > > [hidden email] 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. > > > [[alternative HTML version deleted]] ______________________________________________ [hidden email] 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. |
Free forum by Nabble - Free Resume Builder | Edit this page |