[R] Minimizing a Function with three Parameters
Uwe Ligges
ligges at statistik.uni-dortmund.de
Fri Dec 2 09:04:19 CET 2005
voodooochild at gmx.de wrote:
> 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.
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
> 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
More information about the R-help
mailing list