[R] loop in optim

EdBo n.bowora at gmail.com
Wed Jul 6 04:07:06 CEST 2011


Hi Josh

I have run the code and the structure of the output is what I wanted.
However, the code is giving an identical result for all runs.

I have attached the code I ran below as well as the output. I have just
changed number of runs to match with the size of the data.

a=read.table("D:/hope.txt",header=T)
attach(a)
a
 #likilihood 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) 
#looping now
runs=133/20+1 #total data points divided by number od days in each quater+1
out <- matrix(NA, nrow = runs, ncol = 4,
 dimnames = list(paste("Quater:", 1:runs, sep = ''),
 c("al_j", "au_j", "sigma_j", "b_j"))) 

 for (i in 1:runs) { 
   a[seq(20 * (i - 1) +1, 20 * i), ] 
   out[i, ] <- optim(llik, par = start.par, method = "Nelder-Mead")[[1]] 
 } 
out
#results I am getting
> out
               al_j       au_j       sigma_j      b_j
Quater:1 0.04001525 0.06006251 -7.171336e-25 1.049982
Quater:2 0.04001525 0.06006251 -7.171336e-25 1.049982
Quater:3 0.04001525 0.06006251 -7.171336e-25 1.049982
Quater:4 0.04001525 0.06006251 -7.171336e-25 1.049982
Quater:5 0.04001525 0.06006251 -7.171336e-25 1.049982
Quater:6 0.04001525 0.06006251 -7.171336e-25 1.049982
Quater:7 0.04001525 0.06006251 -7.171336e-25 1.049982
> 







--
View this message in context: http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3647549.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list