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: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of Gene Leynes
Sent: Tuesday, April 07, 2009 12:17 PM
To: r-help at 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 at 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.