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 Thank you in advance. Edward UCT -- View this message in context: http://r.789695.n4.nabble.com/Hint-improve-my-code-tp3641354p3641354.html Sent from the R help mailing list archive at Nabble.com.
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 -- View this message in context: http://r.789695.n4.nabble.com/Hint-improve-my-code-tp3641354p3641520.html Sent from the R help mailing list archive at Nabble.com.
Hi Edward,
One hint to improve your code is simply stylistic. Everyone tends to
have their own favorite way, but adding spaces and line breaks can
help make code much easier to read. Also, you seemed to have used
brackets "[" instead of parentheses "(" for the first two
calls to
sum(). As Berend suggested, some more clarification on what is
supposed to be happening will help you get more feedback. Lastly,
although your data is readable as is, a very convenient way to provide
the list with data is using the function, dput(). Here is a little
example demonstrating its use:
> x <- data.frame(x1 = c("c", "b", "a"),
+ x2 = factor(c("c", "b", "a")),
stringsAsFactors = FALSE)> x ## print to console, two columns appear identical
x1 x2
1 c c
2 b b
3 a a> ## but looking at the str()ucture, they are different
> str(x)
'data.frame': 3 obs. of 2 variables:
$ x1: chr "c" "b" "a"
$ x2: Factor w/ 3 levels "a","b","c": 3 2
1> ## use dput() to provide copy-and-pastable output
> ## that others can use to access the same data as you easy peasy
> dput(x)
structure(list(x1 = c("c", "b", "a"), x2 =
structure(c(3L, 2L,
1L), .Label = c("a", "b", "c"), class =
"factor")), .Names = c("x1",
"x2"), row.names = c(NA, -3L), class =
"data.frame")>
> ## now someone else is ready to go with your data
> newx <- structure(list(x1 = c("c", "b",
"a"), x2 = structure(c(3L, 2L,
+ 1L), .Label = c("a", "b", "c"), class =
"factor")), .Names = c("x1",
+ "x2"), row.names = c(NA, -3L), class =
"data.frame")>
> str(newx)
'data.frame': 3 obs. of 2 variables:
$ x1: chr "c" "b" "a"
$ x2: Factor w/ 3 levels "a","b","c": 3 2
1>
And here is an example of how I might have written your function
(stylistically speaking). I am not advocating this as the "best",
"correct", or "ideal", way. Just giving another example
that (I at
least) find easier to read and make sense of.
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)))
}
}
Cheers,
Josh
On Sat, Jul 2, 2011 at 5:52 PM, EdBo <n.bowora at gmail.com>
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
>
> Thank you in advance.
>
> Edward
> UCT
>
> --
> View this message in context:
http://r.789695.n4.nabble.com/Hint-improve-my-code-tp3641354p3641354.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
http://www.joshuawiley.com/
Please also consider the code below. My attempt to loop through the rows of R_j so that I do not force if() to work on a vector.> llik <- function(par, R_j, R_m) {+ al_j <- par[1] + au_j <- par[2] + sigma_j <- par[3] + b_j <- par[4] + + n=2 + runs=5 + est1=matrix(0,nrow=runs) + start.par=c(al_j=0,au_j=0,sigma_j=0.01,b_j=1) + out1=optim(par=start.par,llik, R_j=R_j, R_m=R_m) + for (i in 1: runs) + { + index_start=2*(i-1)+1 + index_end= 2*i + if(R_j[i]< 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[i]>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[i]==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))) + } + est1[i] <- out1[index_start:index_end] + } + }> est1Error: object 'est1' not found Thank you in advance -- View this message in context: http://r.789695.n4.nabble.com/Hint-improve-my-code-tp3641354p3642707.html Sent from the R help mailing list archive at Nabble.com.