maximizing a function

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

maximizing a function

Nigel.Walker
Hi list,
I would like to maximize a likelihood function to obtain estimates. The problem is that i first have to have functions other functions that are in the likelihood that depend on the parameters. I was thinking of using optim but i do not know how to go about starting. A example of what i wish to do is given below. The parameters i wish to estimate are phi and gamma.

How would i go about maximizing the likelihood function to obtain estimates?

Any help would be much appreciated

Nigel

# initial starting points
phi<-0.2
gamma<-0.1
n<-200
h<-c()
h[1]<-1
e<-c()
e[1]<-1
for (i in 2:n){
h[i]<-phi*e[i-1]+gamma*h[i-1]+3
e[i]<-rnorm(1,0,h[i])
}
likelihood<- sum(log(h)+e^2/2*h)

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: maximizing a function

Uwe Ligges


[hidden email] wrote:

> Hi list,
> I would like to maximize a likelihood function to obtain estimates. The problem is that i first have to have functions other functions that are in the likelihood that depend on the parameters. I was thinking of using optim but i do not know how to go about starting. A example of what i wish to do is given below. The parameters i wish to estimate are phi and gamma.
>
> How would i go about maximizing the likelihood function to obtain estimates?
>
> Any help would be much appreciated
>
> Nigel
>
> # initial starting points
> phi<-0.2
> gamma<-0.1
> n<-200
> h<-c()
> h[1]<-1
> e<-c()
> e[1]<-1
> for (i in 2:n){
> h[i]<-phi*e[i-1]+gamma*h[i-1]+3
> e[i]<-rnorm(1,0,h[i])
> }
> likelihood<- sum(log(h)+e^2/2*h)



Please read ?optim more precisely.

You want

lik <- function(par, n=200){
     phi <- par[1]
     gamma <- par[2]
     h <- numeric(n)
     h[1] <- 1
     e <- numeric(n)
     e[1] <- 1
     for (i in 2:n){
         h[i] <- phi*e[i-1] + gamma*h[i-1] + 3
         e[i] <- rnorm(1,0,h[i])
     }
     ## optim minimizes, hence using negative sum:
     -sum(log(h) + e^2/2*h)
}

set.seed(1)
optim(c(phi=0.3, gamma=0.1), lik, n=200)


and you will find that your function is highly dependend on the random
numbers you are generating, hence using a seed. I wonder where the
random term comes from. Probably there is some error in the likelihood?

Uwe Ligges

______________________________________________
[hidden email] 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.