Very nice trick!!! Thank you!
When you combine your trick with the "nls" function, you get the whole
thing:
#(I have also simplified the model since the first post)
###################################################################
#Original Method: uses optim, but doesn't create "model":
###################################################################
x=matrix(c(0.04, 0.08, 0.01, 0.0398, 0.081, 0.00915, 0.04057, 0.07944,
0.00994, 0.04137, 0.07949, 0.01132, 0.0421, 0.08273, 0.00947, 0.04237,
0.08058, 0.00969, 0.0425, 0.07952, 0.00919, 0.04305, 0.07717, 0.00908,
0.04319, 0.07576, 0.0061, 0.04298, 0.07557, 0.00539, 0.04287, 0.07244,
0.00542, 0.04318, 0.071, 0.00388, 0.04348, 0.07053, 0.00375, 0.04335,
0.07075, 0.00364, 0.04386, 0.07481, 0.00296),ncol=3,byrow=T)
x=cbind(1,x)
y=matrix(c(0.02847,0.02831,0.03255,0.02736,0.03289,0.02774,0.03192,0.03141,
0.02927,0.02648,0.02807,0.03047,0.03046,0.02541,0.02729))
plot(x=1:15,ylim=c(max(x[,-1]),min(x[,-1])))
matlines(cbind(x[,-1],y))
# "normalization" fun to make the rest of the x's add up to 1
normalize=function(x)c(x[1], x[2:length(x)]/sum(x[2:length(x)]))
# objective function:
fn=function(ws){
ws=normalize(ws)
sum((x %*% ws - yhat)^2)}
llim = c(-Inf,rep(0,ncol(x)-1)) # alpha is -Inf to Inf
ulim = c( Inf,rep(1,ncol(x)-1)) # betas are 0 to 1
th=matrix(c(0,rep(1/3,3)))
o = optim(th, fn, method="L-BFGS-B",lower=llim, upper=ulim)
NormPar = normalize(o$par)
cat('pars:', NormPar,'\n', 'tot:', sum(NormPar), '\n')
###################################################################
# Dr. Varadhan's Method to make betas sum to 1:
###################################################################
xnew=rbind(c(0, 1, 1, 1),x)
ynew=rbind(1,y)
ans <- lm(ynew ~ xnew - 1)
summary(ans)
sum(ans$resid^2) # compare with the objective function value from
cat('pars:', ans$coef,'\n', 'tot:', sum(ans$coef[2:4]), '\n')
###################################################################
# Dr. Varadhan's Method WITH nls Method
###################################################################
xnew=rbind(c(0, 1, 1, 1),x)
ynew=rbind(1,y)
df=data.frame(cbind(ynew,xnew))
colnames(df)=c('y','x0','x1','x2','x3')
aformula = y ~ b0*x0 + b1*x1 + b2*x2 + b3*x3
ans2 <- nls(aformula,data=df,,algorithm='port',
start=list(b0=0,b1=1/3,b2=1/3,b3=1/3),
lower=c(-Inf,rep(0,3)),upper=c(Inf,rep(1,3)),trace=T)
summary(ans2)
sum(residuals(ans2)^2) # compare with the objective function value from
cat('pars:', coef(ans2),'\n', 'tot:', sum(coef(ans2)[2:4]), '\n')
On Tue, Apr 7, 2009 at 12:37 PM, Ravi Varadhan <rvarad...@jhmi.edu> wrote:
> Hi Gene,
>
> Try the following approach using lm():
>
> set.seed(121)
> x1=.04
> for (i in 1:14) x1[i+1]=x1[i]*(1+rnorm(1)*.008)+.00025
> x2=.08
> for (i in 1:14) x2[i+1]=x2[i]*(1+rnorm(1)*.03)-.0018
> x3=.01
> for (i in 1:14) x3[i+1]=x3[i]*(1+rnorm(1)*.15)-.0008
>
> b=matrix(c(0.6,0.0,0.4)) # why don't you assume a value for intercept?
> x=matrix(cbind(x1,x2,x3),ncol=3)
>
> y=x%*%b # the 'real' y
> yhat=y+runif(15)*.006 # the observed y
> x=cbind(1,x)
>
> # here is the simple trick: add a row to X matrix and an element to yhat
> vector to impose constraints
> #
> xadd <- c(0, 1, 1, 1)
> xnew <- rbind(xadd, x)
> ynew <- c(1, yhat)
>
> ans <- lm(ynew ~ xnew - 1)
> summary(ans)
> sum(ans$coef[2:4])
> sum(ans$resid^2) # compare with the objective function value from
> your approach
> > sum(ans$resid^2)
> [1] 2.677474e-05
> >
>
> > o$value
> [1] 2.718646e-05
> >
>
> Note that the sum of squared residuals from lm() is smaller than your value
> from optim(). Although, this approach works well in your example, it does
> not guarantee that the coefficients are between 0 and 1.
>
> Ravi.
>
>
> ----------------------------------------------------------------------------
> -------
>
> Ravi Varadhan, Ph.D.
>
> Assistant Professor, The Center on Aging and Health
>
> Division of Geriatric Medicine and Gerontology
>
> Johns Hopkins University
>
> Ph: (410) 502-2619
>
> Fax: (410) 614-9625
>
> Email: rvarad...@jhmi.edu
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
>
> ----------------------------------------------------------------------------
> --------
>
>
> -----Original Message-----
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> On
> Behalf Of Gene Leynes
> Sent: Tuesday, April 07, 2009 12:17 PM
> To: r-help@r-project.org
> Subject: [R] get optim results into a model object
>
> Hello all, I have an optimization routine that is giving me good results,
> but the results are not in the nice "model" format like "lm". How can I
> get
> optim results into a model so that I can use the clever 'fitted',
> 'residuals', and 'summary' functions?
>
> Using optim is the only way that I was able to make a model that
> 1) sums the betas to 1,
> 2) constrains the betas to positive numbers less than 1
> 3) does not constrain alpha
> (The constrOptim cousin wasn't very accurate, and was very slow.)
>
> Here is an example of some code, the results of which I would like to get
> into a model with the form y ~ alpha + REALPAR * x where 'REALPAR' is the
> "normalized" output at the very end
>
> many thanks!!!!
> ++++++++++++++++Code Below++++++++++++++++++++++++++++
>
> set.seed(121)
> x1=.04
> for (i in 1:14) x1[i+1]=x1[i]*(1+rnorm(1)*.008)+.00025
> x2=.08
> for (i in 1:14) x2[i+1]=x2[i]*(1+rnorm(1)*.03)-.0018
> x3=.01
> for (i in 1:14) x3[i+1]=x3[i]*(1+rnorm(1)*.15)-.0008
>
> b=matrix(c(0.6,0.0,0.4))
> x=matrix(cbind(x1,x2,x3),ncol=3)
> y=x%*%b # the 'real' y
> yhat=y+runif(15)*.006 # the observed y
>
> plot(x=1:15,ylim=c(min(x1,x2,x3),max(x1,x2,x3)))
> matlines(cbind(x,y,yhat))
>
> # Add a constant to x (for alpha)
> x=cbind(1,x)
> # "normalization" fun to make the rest of the x's add up to 1
> normalize=function(x)c(x[1], x[2:length(x)]/sum(x[2:length(x)]))
> # objective function:
> fn=function(ws){
> ws=normalize(ws)
> sum((x %*% ws - yhat)^2)}
>
> llim = c(-Inf,rep(0,ncol(x)-1)) # alpha (col 1 of 'x') is -Inf to Inf ulim
> = c( Inf,rep(1,ncol(x)-1)) # betas (cols 2:4 of 'x') are 0 to 1
> th=matrix(c(0,rep(1/3,3)))
> o = optim(th, fn, method="L-BFGS-B",lower=llim, upper=ulim) o REALPAR =
> normalize(o$par) REALPAR
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> 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.
>
>
[[alternative HTML version deleted]]
______________________________________________
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.
