[R] loop in optim

Joshua Wiley jwiley.psych at gmail.com
Thu Jul 7 07:00:58 CEST 2011


I don't know what else to say.  Your code looks right to me, and it
all runs.  I would check the value of a at each loop:

for (i in 1:4) {
  a <- afull[seq(20 * (i - 1) +1, 20 * i), ]
  print(a) # so you can see what it is
  out[i, ] <- optim(llik, par = start.par, method = "Nelder-Mead")[[1]]
}

Also, I would try running the code in a clean version of R.  That is,
a version of R without any non-standard packages loaded, or clutter in
the workspace.  You can do this by shutting down R, deleting the old
workspace, and then starting a new session (there are many other ways
to do the same thing, that's just one).

Josh

On Wed, Jul 6, 2011 at 4:34 PM, EdBo <n.bowora at gmail.com> wrote:
> I am sorry if I sound stupid but I am not able to correct the error
> even after running this code.
>
>> afull=read.table("D:/hope.txt",header=T)
>> 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]]
> + }
>> out
>           al_j      au_j       sigma_j      b_j
> Run:1 0.1088116 0.1621605 -1.554167e-24 0.969153
> Run:2 0.1088116 0.1621605 -1.554167e-24 0.969153
> Run:3 0.1088116 0.1621605 -1.554167e-24 0.969153
> Run:4 0.1088116 0.1621605 -1.554167e-24 0.959875
>
> On 6 July 2011 11:46, Berend Hasselman [via R]
> <ml-node+3648171-330506108-247474 at n4.nabble.com> wrote:
>> EdBo wrote:
>> You are right Joshua.
>>
>> I changed the code because I failed to understand how you attached the full
>> data set. How you made the data part of your code.
>>
>> I am new to R so I am used to one way of attaching data(the way I redone
>> it).
>>
>> You don't need to "attach" the data by using attach().
>> You read the data into an object afull and then select the part you need and
>> store that in object a.
>>
>> BTW: shouldn't the for (i in 1:4) be for (i in 1:3) if I understand the
>> original question correctly?
>>
>> Berend
>>
>> ________________________________
>> If you reply to this email, your message will be added to the discussion
>> below:
>> http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3648171.html
>> To unsubscribe from loop in optim, click here.
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3650297.html
> Sent from the R help mailing list archive at Nabble.com.
>        [[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.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
https://joshuawiley.com/



More information about the R-help mailing list