Optimisation with Normalisation Constraint

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Optimisation with Normalisation Constraint

Lorenzo Isella
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.
Reply | Threaded
Open this post in threaded view
|

Re: Optimisation with Normalisation Constraint

Jeff Newmiller
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.
Reply | Threaded
Open this post in threaded view
|

Re: Optimisation with Normalisation Constraint

David Winsemius
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.
Reply | Threaded
Open this post in threaded view
|

Re: Optimisation with Normalisation Constraint

Hans W Borchers-2
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.