[R] loop in optim
Joshua Wiley
jwiley.psych at gmail.com
Tue Jul 5 19:17:49 CEST 2011
Do you want something like this?
## Your data
afull <- read.table(textConnection("
R_j R_m
-0.0625 0.002320654
0 -0.004642807
0.033333333 0.005936332
0.032258065 0.001060848
0 0.007114057
0.015625 0.005581558
0 0.002974794
0.015384615 0.004215271
0.060606061 0.005073116
0.028571429 -0.006001279
0 -0.002789594
0.013888889 0.00770633
0 0.000371663
0.02739726 -0.004224228
-0.04 0.008362539
0 -0.010951605
0 0.004682924
0.013888889 0.011839993
-0.01369863 0.004210383
-0.027777778 -0.04658949
0 0.00987272
-0.057142857 -0.062203157
-0.03030303 -0.119177639
0.09375 0.077054642
0 -0.022763619
-0.057142857 0.050408775
0 0.024706076
-0.03030303 0.004043701
0.0625 0.004951088
0 -0.005968731
0 -0.038292548
0 0.013381097
0.014705882 0.006424728
-0.014492754 -0.020115626
0 -0.004837891
-0.029411765 -0.022054654
0.03030303 0.008936428
0.044117647 8.16925E-05
0 -0.004827246
-0.042253521 0.004653096
-0.014705882 -0.004222151
0.029850746 0.000107267
-0.028985507 -0.001783206
0.029850746 -0.006372981
0.014492754 0.005492374
-0.028571429 -0.009005846
0 0.001031683
0.044117647 0.002800551"), header = TRUE)
closeAllConnections()
## likelihood function
llik = function(x)
{
al_j=x[1]; au_j=x[2]; sigma_j=x[3]; b_j=x[4]
sum(na.rm=T,
ifelse(a$R_j< 0, log(1/(2*pi*(sigma_j^2)))-
(1/(2*(sigma_j^2))*(a$R_j+al_j-b_j*a$R_m))^2,
ifelse(a$R_j>0 , log(1/(2*pi*(sigma_j^2)))-
(1/(2*(sigma_j^2))*(a$R_j+au_j-b_j*a$R_m))^2,
log(ifelse (( pnorm (au_j, mean=b_j *
a$R_m, sd= sqrt(sigma_j^2))-
pnorm(al_j, mean=b_j * a$R_m, sd=sqrt (sigma_j^2)
)) > 0,
(pnorm (au_j,mean=b_j * a$R_m,
sd= sqrt(sigma_j^2))-
pnorm(al_j, mean=b_j * a$R_m, sd= sqrt(sigma_j^2)
)),
1)) ))
)
}
start.par = c(-0.01,0.01,0.1,1)
out <- matrix(NA, nrow = 4, ncol = 4, dimnames = list(
paste("Run:", 1:4, sep = ''),
c("al_j", "au_j", "sigma_j", "b_j")))
## Estimate parameters based on rows 0-20, 21-40, 41-60 of 'afull'
for (i in 1:4) {
a <- afull[seq(20 * (i - 1) +1, 20 * i), ]
out[i, ] <- optim(llik, par = start.par, method = "Nelder-Mead")[[1]]
}
## Yields
> out
al_j au_j sigma_j b_j
Run:1 0.04001776 0.06010743 1.092618e-24 1.049971
Run:2 0.04002135 0.06008513 -7.156966e-25 1.049976
Run:3 0.04714390 0.27258724 3.303320e-24 0.948988
Run:4 -0.01000000 0.01000000 1.000000e-01 1.000000
On Mon, Jul 4, 2011 at 8:21 PM, EdBo <n.bowora at gmail.com> wrote:
> Hi
>
> I have re-worked on my likelihood function and it is now working(#the code
> is below#).
>
> May you help me correct my loop function.
>
> I want optim to estimates al_j; au_j; sigma_j; b_j by looking at 0 to 20,
> 21 to 40, 41 to 60 data points.
>
> The final result should have 4 columns of each of the estimates AND 4 rows
> of each of 0 to 20, 21 to 40, 41 to 60.
>
> #likelihood function
> a=read.table("D:/hope.txt",header=T)
> attach(a)
> a
> llik = function(x)
> {
> al_j=x[1]; au_j=x[2]; sigma_j=x[3]; b_j=x[4]
> sum(na.rm=T,
> ifelse(a$R_j< 0, log(1/(2*pi*(sigma_j^2)))-
> (1/(2*(sigma_j^2))*(a$R_j+al_j-b_j*a$R_m))^2,
> ifelse(a$R_j>0 , log(1/(2*pi*(sigma_j^2)))-
> (1/(2*(sigma_j^2))*(a$R_j+au_j-b_j*a$R_m))^2,
> log(ifelse (( pnorm (au_j, mean=b_j * a$R_m, sd= sqrt(sigma_j^2))-
> pnorm(al_j, mean=b_j * a$R_m, sd=sqrt (sigma_j^2)
> )) > 0,
> (pnorm (au_j,mean=b_j * a$R_m, sd= sqrt(sigma_j^2))-
> pnorm(al_j, mean=b_j * a$R_m, sd= sqrt(sigma_j^2)
> )),
> 1)) ))
> )
> }
> start.par = c(-0.01,0.01,0.1,1)
> out1 = optim(llik, par=start.par, method="Nelder-Mead")
> out1
>
> --
> View this message in context: http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3645031.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> 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.
>
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
https://joshuawiley.com/
More information about the R-help
mailing list