Thank you for your kind responses.
What I want my likelihood function to do is to look at each row and
see if R_j is (> 0 or < 0 or == 0) then decide which of the three
parts of the likelihood function apply to that row.
At the end, it should use all that information to calculate the best
estimate of each of the parameters.
On 3 July 2011 06:59, Berend Hasselman [via R]
<ml-node+3641520-2044208009-247...@n4.nabble.com> wrote:
> EdBo wrote:
> Hi
>
> I have developed the code below. I am worried that the parameters I want to
> be estimated are "not being found" when I ran my code. Is there a way I can
> code them so that R recognize that they should be estimated.
>
> This is the error I am getting.
>
>> out1=optim(llik,par=start.par)
> Error in pnorm(au_j, mean = b_j * R_m, sd = sigma_j) :
> Â object 'au_j' not found
>
> #Yet al_j,au_j,sigma_j and b_j are just estimates that balance the
> likelihood function?
>
> llik=function(R_j,R_m)
> if(R_j< 0)
> {
> sum[log(1/(2*pi*(sigma_j^2)))-(1/(2*(sigma_j^2))*(R_j+al_j-b_j*R_m))^2]
> }else if(R_j>0)
> {
> sum[log(1/(2*pi*(sigma_j^2)))-(1/(2*(sigma_j^2))*(R_j+au_j-b_j*R_m))^2]
> }else if(R_j==0)
> {
> sum(log(pnorm(au_j,mean=b_j*R_m,sd=sigma_j)-pnorm(al_j,mean=b_j*R_m,sd=sigma_j)))
> }
> start.par=c(al_j=0,au_j=0,sigma_j=0.01,b_j=1)
> out1=optim(par=start.par,llik)
>
> My Data
>
>    R_j     R_m
> Â 2e-03 Â 0.026567295
> Â 3e-03 Â 0.009798475
> Â 5e-02 Â 0.008497274
> Â -1e-02 Â 0.012464578
> Â -9e-04 Â 0.002896023
> Â 9e-02 Â 0.000879473
> Â 1e-02 Â 0.003194435
> Â 6e-04 Â 0.010281122
>
> You are only passing the R_j and R_m data as argument to your function.
> The parameters that are to be used for the maximization must also be passed
> as arguments.
> This is clearly explained in the help page for optim. So your function
> should read
>
> llik <- function(par, R_j, R_m) {
>   al_j    <- par[1]
>   au_j    <- par[2]
> Â Â sigma_j <- par[3]
>   b_j    <- par[4]
>
> Â Â if(R_j< 0) {
>
> sum(log(1/(2*pi*(sigma_j^2)))-(1/(2*(sigma_j^2))*(R_j+al_j-b_j*R_m))^2)
> Â Â } else if(R_j>0) Â {
>
> sum(log(1/(2*pi*(sigma_j^2)))-(1/(2*(sigma_j^2))*(R_j+au_j-b_j*R_m))^2)
> Â Â } else if(R_j==0) {
>
> sum(log(pnorm(au_j,mean=b_j*R_m,sd=sigma_j)-pnorm(al_j,mean=b_j*R_m,sd=sigma_j)))
> Â Â }
> }
>
> And then use optim as follows:
>
> out1 <- optim(par=start.par,llik, R_j=R_j, R_m=R_m)
>
>
> But this is not going to work as you want it to.
> You'll get warnings:
>
> 1: In if (R_j < 0) { ... :
> Â the condition has length > 1 and only the first element will be used
>
> R_j is a vector and if() takes a scalar (or vector of length 1).
>
> So it is not at all clear what the purpose of the if statement is.
>
> Berend
>
>
>
> ________________________________
> If you reply to this email, your message will be added to the discussion
> below:
> http://r.789695.n4.nabble.com/Hint-improve-my-code-tp3641354p3641520.html
> To unsubscribe from Hint improve my code, click here.
--
View this message in context:
http://r.789695.n4.nabble.com/Hint-improve-my-code-tp3641354p3642580.html
Sent from the R help mailing list archive at Nabble.com.
[[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.
