R help - Dec 2005 - Minimizing a Function with three Parameters

Hi,

I'm trying to get maximum likelihood estimates of \alpha, \beta_0 and 
\beta_1, this can be achieved by solving the following three equations:

n / \alpha + \sum\limits_{i=1}^{n} ln(\psihat(i)) - 
\sum\limits_{i=1}^{n} ( ln(x_i + \psihat(i)) ) = 0

\alpha \sum\limits_{i=1}^{n} 1/(psihat(i)) - (\alpha+1) 
\sum\limits_{i=1}^{n} ( 1 / (x_i + \psihat(i)) ) = 0

\alpha \sum\limits_{i=1}^{n} ( i / \psihat(i) ) - (\alpha + 1) 
\sum\limits_{i=1}^{n} ( i / (x_i + \psihat(i)) ) = 0

where \psihat=\beta_0 + \beta_1 * i. Now i want to get iterated values 
for \alpha, \beta_0 and \beta_1, so i used the following implementation

# first equation
l1 <- function(beta0,beta1,alpha,x) {
  n<-length(x)
  s2<-length(x)
    for(i in 1:n) {
    s2[i]<-log(beta0+beta1*i)-log(x[i]+beta0+beta1*i)
    }
  s2<-sum(s2)
  return((n/alpha)+s2)
}

# second equation
l2 <- function(beta0,beta1,alpha,x) {
  n<-length(x)
  s1<-length(x)
  s2<-length(x)
    for(i in 1:n) {
    s1[i]<-1/(beta0+beta1*i)
    s2[i]<-1/(beta0+beta1*i+x[i])
    }
  s1<-sum(s1)
  s2<-sum(s2)
  return(alpha*s1-(alpha+1)*s2)
}

#third equation
l3 <- function(beta0,beta1,alpha,x) {
  n<-length(x)
  s1<-length(x)
  s2<-length(x)
    for(i in 1:n) {
    s1[i]<-i/(beta0+beta1*i)
    s2[i]<-i/(x[i]+beta0+beta1*i)
    }
  s1<-sum(s1)
  s2<-sum(s2)
  return(alpha*s1-(alpha+1)*s2)
}

# all equations in one
gl <- function(beta0,beta1,alpha,x) {
  l1(beta0,beta1,alpha,x)^2 + l2(beta0,beta1,alpha,x)^2 + 
l3(beta0,beta1,alpha,x)^2
}

#iteration with optim
optim(c(1,1,1),gl,x)

i get always an error massage. Is optim anyway the 'right' method to get
all three parameters iterated at the same time?

best regards
Andreas

voodooochild at gmx.de wrote:
> Hi,
> 
> I'm trying to get maximum likelihood estimates of \alpha, \beta_0 and 
> \beta_1, this can be achieved by solving the following three equations:
> 
> n / \alpha + \sum\limits_{i=1}^{n} ln(\psihat(i)) - 
> \sum\limits_{i=1}^{n} ( ln(x_i + \psihat(i)) ) = 0
> 
> \alpha \sum\limits_{i=1}^{n} 1/(psihat(i)) - (\alpha+1) 
> \sum\limits_{i=1}^{n} ( 1 / (x_i + \psihat(i)) ) = 0
> 
> \alpha \sum\limits_{i=1}^{n} ( i / \psihat(i) ) - (\alpha + 1) 
> \sum\limits_{i=1}^{n} ( i / (x_i + \psihat(i)) ) = 0
> 
> where \psihat=\beta_0 + \beta_1 * i. Now i want to get iterated values 
> for \alpha, \beta_0 and \beta_1, so i used the following implementation
> 
> # first equation
> l1 <- function(beta0,beta1,alpha,x) {
>   n<-length(x)
>   s2<-length(x)
>     for(i in 1:n) {
>     s2[i]<-log(beta0+beta1*i)-log(x[i]+beta0+beta1*i)
>     }
>   s2<-sum(s2)
>   return((n/alpha)+s2)
> }
> 
> # second equation
> l2 <- function(beta0,beta1,alpha,x) {
>   n<-length(x)
>   s1<-length(x)
>   s2<-length(x)
>     for(i in 1:n) {
>     s1[i]<-1/(beta0+beta1*i)
>     s2[i]<-1/(beta0+beta1*i+x[i])
>     }
>   s1<-sum(s1)
>   s2<-sum(s2)
>   return(alpha*s1-(alpha+1)*s2)
> }
> 
> #third equation
> l3 <- function(beta0,beta1,alpha,x) {
>   n<-length(x)
>   s1<-length(x)
>   s2<-length(x)
>     for(i in 1:n) {
>     s1[i]<-i/(beta0+beta1*i)
>     s2[i]<-i/(x[i]+beta0+beta1*i)
>     }
>   s1<-sum(s1)
>   s2<-sum(s2)
>   return(alpha*s1-(alpha+1)*s2)
> }
> 
> # all equations in one
> gl <- function(beta0,beta1,alpha,x) {
>   l1(beta0,beta1,alpha,x)^2 + l2(beta0,beta1,alpha,x)^2 + 
> l3(beta0,beta1,alpha,x)^2
> }
> 
> #iteration with optim
> optim(c(1,1,1),gl,x)
> 
> i get always an error massage. Is optim anyway the 'right' method
to get
> all three parameters iterated at the same time?
> 
> best regards
> Andreas
> 
Hi, Andreas,

You've misread the help file for ?optim.

  fn: A function to be minimized (or maximized), with first
           argument the vector of parameters over which minimization is
           to take place.  It should return a scalar result.

So, your objective function should look like

gl <- function(beta, x) {
   beta0 <- beta[1]
   beta1 <- beta[2]
   alpha <- beta[3]
   v1 <- l1(beta0, beta1, alpha, x)^2
   v2 <- l2(beta0, beta1, alpha, x)^2
   v3 <- l3(beta0, beta1, alpha, x)^2
   v1 + v2 + v3
}

Also, are you aware of ?mle in the stats4 package?

HTH,

--sundar

voodooochild at gmx.de wrote:
>Hi,
>
>I'm trying to get maximum likelihood estimates of \alpha, \beta_0 and 
>\beta_1, this can be achieved by solving the following three equations:
>
>n / \alpha + \sum\limits_{i=1}^{n} ln(\psihat(i)) - 
>\sum\limits_{i=1}^{n} ( ln(x_i + \psihat(i)) ) = 0
>
>\alpha \sum\limits_{i=1}^{n} 1/(psihat(i)) - (\alpha+1) 
>\sum\limits_{i=1}^{n} ( 1 / (x_i + \psihat(i)) ) = 0
>
>\alpha \sum\limits_{i=1}^{n} ( i / \psihat(i) ) - (\alpha + 1) 
>\sum\limits_{i=1}^{n} ( i / (x_i + \psihat(i)) ) = 0
>
>where \psihat=\beta_0 + \beta_1 * i. Now i want to get iterated values 
>for \alpha, \beta_0 and \beta_1, so i used the following implementation
>
># first equation
>l1 <- function(beta0,beta1,alpha,x) {
>  n<-length(x)
>  s2<-length(x)
>    for(i in 1:n) {
>    s2[i]<-log(beta0+beta1*i)-log(x[i]+beta0+beta1*i)
>    }
>  s2<-sum(s2)
>  return((n/alpha)+s2)
>}
>
># second equation
>l2 <- function(beta0,beta1,alpha,x) {
>  n<-length(x)
>  s1<-length(x)
>  s2<-length(x)
>    for(i in 1:n) {
>    s1[i]<-1/(beta0+beta1*i)
>    s2[i]<-1/(beta0+beta1*i+x[i])
>    }
>  s1<-sum(s1)
>  s2<-sum(s2)
>  return(alpha*s1-(alpha+1)*s2)
>}
>
>#third equation
>l3 <- function(beta0,beta1,alpha,x) {
>  n<-length(x)
>  s1<-length(x)
>  s2<-length(x)
>    for(i in 1:n) {
>    s1[i]<-i/(beta0+beta1*i)
>    s2[i]<-i/(x[i]+beta0+beta1*i)
>    }
>  s1<-sum(s1)
>  s2<-sum(s2)
>  return(alpha*s1-(alpha+1)*s2)
>}
>
># all equations in one
>gl <- function(beta0,beta1,alpha,x) {
>  l1(beta0,beta1,alpha,x)^2 + l2(beta0,beta1,alpha,x)^2 + 
>l3(beta0,beta1,alpha,x)^2
>}
>
>#iteration with optim
>optim(c(1,1,1),gl,x)
>
>i get always an error massage. Is optim anyway the 'right' method to
get
>all three parameters iterated at the same time?
>
>best regards
>Andreas
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>
>
>  
>
hi sundar,

your advice has helped very much, thanks a lot.

now i have another model where instead of i    i2 is used, but i don't 
now way i got so large estimates?

x<-c(10,8,14,17,15,22,19,27,35,40)

# 1.Gleichung

l1 <- function(beta0,beta1,alpha,x) {
 n<-length(x)
 s2<-length(x)
   for(i in 1:n) {
   s2[i]<-log(beta0+beta1*i2)-log(x[i]+beta0+beta1*i2)
   }
 s2<-sum(s2)
 return((n/alpha)+s2)
}


# 2.Gleichung

l2 <- function(beta0,beta1,alpha,x) {
 n<-length(x)
 s1<-length(x)
 s2<-length(x)
   for(i in 1:n) {
   s1[i]<-1/(beta0+beta1*i2)
   s2[i]<-1/(beta0+beta1*i2+x[i])
   }
 s1<-sum(s1)
 s2<-sum(s2)
 return(alpha*s1-(alpha+1)*s2)
}

# 3.Gleichung

l3 <- function(beta0,beta1,alpha,x) {
 n<-length(x)
 s1<-length(x)
 s2<-length(x)
   for(i in 1:n) {
   s1[i]<-(i2)/(beta0+beta1*i2)
   s2[i]<-(i2)/(x[i]+beta0+beta1*i2)
   }
 s1<-sum(s1)
 s2<-sum(s2)
 return(alpha*s1-(alpha+1)*s2)
}

# Zusammenf??gen aller Teile

gl <- function(beta,x) {
 beta0<-beta[1]
 beta1<-beta[2]
 alpha<-beta[3]
 v1<-l1(beta0,beta1,alpha,x)2
 v2<-l2(beta0,beta1,alpha,x)2
 v3<-l3(beta0,beta1,alpha,x)2
 v1+v2+v3
}


# Nullstellensuche mit Nelder-Mead

optim(c(20000,6000,20000),gl,x=x,control=list(reltol=1e-12))

the values should be alpha=20485, beta0=19209 and beta1=6011

and another point is, what is a good method to find good starting values 
for 'optim'. it seems, that i only get the desired values when the 
starting values are in the same region. I used 
control=list(reltol=1e-12), but it seems, that then it is also important 
to have the starting values in the same region as the the desired values.

regards
andreas