Skip to content

Minimizing a Function with three Parameters

4 messages · voodooochild@gmx.de, Sundar Dorai-Raj, Uwe Ligges

#
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, 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 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
#
voodooochild at gmx.de wrote:

            
Depends on the response surface. If the latter behaves well, you can be 
lucky, if it does behave really ill, you will get grey hair during your 
further research on this topic ...
Please read some textbook on numerical optimization.

Uwe Ligges