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.
