Dear All,
I have a problem I haver been struggling with for a while: I need to carry out a non-linear fit (and this is the easy part). I have a set of discrete values {x1,x2...xN} and the corresponding {y1, y2...yN}. The difficulty is that I would like the linear fit to preserve the sum of the values y1+y2+...yN. I give an example below (for which there may even be an analytical solution, but that is not the point here) ############################################################################ library(minpack.lm) set.seed(124) z <- rexp(3000,3) zf <- z[z<= 0.5 | z>=0.9] myhist <- hist(zf, plot=FALSE) df <- data.frame(x=myhist$mids, y=myhist$density) myfit <- nlsLM(y~(A*exp(-lambda*x)) ,data=df, start=list(A=1,lambda=1)) > sum(myhist$density) [1] 5 > sum(predict(myfit)) [1] 4.931496 ############################################################################ I would like sum(predict(myfit)) to be exactly 5 from the start, without renormalising a posteriori the fit. Any suggestion is appreciated. Cheers Lorenzo ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see 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. |
I recommend posting this on a mathematics discussion forum like Stack Exchange and (re-)reading the Posting Guide for this mailing list.
I think you are going to need to re-write your model function to algebraically combine your original model along with the constraint, and then use the original model alone for prediction... but I haven't tried it so might be quite far off the mark. On June 20, 2018 8:50:48 AM PDT, Lorenzo Isella <[hidden email]> wrote: >Dear All, >I have a problem I haver been struggling with for a while: I need to >carry out a non-linear fit (and this is the >easy part). >I have a set of discrete values {x1,x2...xN} and the corresponding >{y1, y2...yN}. The difficulty is that I would like the linear fit to >preserve the sum of the values y1+y2+...yN. >I give an example below (for which there may even be an analytical >solution, but that is not the point here) > >############################################################################ >library(minpack.lm) > > > >set.seed(124) > >z <- rexp(3000,3) > > >zf <- z[z<= 0.5 | z>=0.9] > >myhist <- hist(zf, plot=FALSE) > > > >df <- data.frame(x=myhist$mids, y=myhist$density) > > > >myfit <- nlsLM(y~(A*exp(-lambda*x)) > ,data=df, start=list(A=1,lambda=1)) > > > >> sum(myhist$density) >[1] 5 >> sum(predict(myfit)) >[1] 4.931496 > >############################################################################ >I would like sum(predict(myfit)) to be exactly 5 from the start, >without renormalising a posteriori the fit. > >Any suggestion is appreciated. >Cheers > >Lorenzo > >______________________________________________ >[hidden email] mailing list -- To UNSUBSCRIBE and more, see >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. -- Sent from my phone. Please excuse my brevity. ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see 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 Lorenzo Isella
> On Jun 20, 2018, at 8:50 AM, Lorenzo Isella <[hidden email]> wrote: > > Dear All, > I have a problem I haver been struggling with for a while: I need to > carry out a non-linear fit (and this is the > easy part). > I have a set of discrete values {x1,x2...xN} and the corresponding > {y1, y2...yN}. The difficulty is that I would like the linear fit to > preserve the sum of the values y1+y2+...yN. > I give an example below (for which there may even be an analytical > solution, but that is not the point here) > > ############################################################################ > library(minpack.lm) > > > > set.seed(124) > > z <- rexp(3000,3) > > > zf <- z[z<= 0.5 | z>=0.9] > > myhist <- hist(zf, plot=FALSE) > > > df <- data.frame(x=myhist$mids, y=myhist$density) > > > > myfit <- nlsLM(y~(A*exp(-lambda*x)) > ,data=df, start=list(A=1,lambda=1)) > > > >> sum(myhist$density) > [1] 5 >> sum(predict(myfit)) > [1] 4.931496 > > ############################################################################ > I would like sum(predict(myfit)) to be exactly 5 from the start, > without renormalising a posteriori the fit. Wouldn't that happen if you minimized that absolute deviations from the fit rather than minimizing the sums of squares?? > > Any suggestion is appreciated. > Cheers > > Lorenzo > David Winsemius Alameda, CA, USA 'Any technology distinguishable from magic is insufficiently advanced.' -Gehm's Corollary to Clarke's Third Law ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see 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 Lorenzo Isella
One way will be to solve this as an ordinary optimization problem with
an equality constraint. Function `alabama::auglag` can do this: library(alabama) fn <- function(p) sum((df$y - p[1]*exp(-p[2]*df$x))^2) heq <- function(p) sum(p[1]*exp(-p[2]*df$x)) - 5 # Start with initial values near first solution sol <- auglag(c(4, 4), fn=fn, heq=heq) Solution with myfit: 4.134 4.078 Solution with auglag: 4.126763 4.017768 The equality constraint is still fulfilled: a <- sol$par[1]; lambda <- sol$par[2] sum(a*exp(-lambda*df$x)) ## [1] 5 Plot the differences of these two solutions: plot(df$x, a*exp(-lambda*df$x) - predict(myfit), type='l') Interestingly, the difference between 5 and `predict(myfit)` is *not* just evenly spread across all 15 points. ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see 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 | Edit this page |