Hi all,
I apply "optim" to the function "obj", which minimizes the
goodness of
fit statistic and obtains Pearson minimum chi-squared estimate for x[1],
x[2] and x[3]. The vector "fr" contains the four observed frequencies.
Since "fr[i]" appears in the denominator, I would like to substitute
"0"
in the sum if fr[i]=0.
I tried an "ifelse" condition which works, but in some cases the
solution seems rather odd. Is there anything wrong with the way second
objective function is coded?
obj<-function(x){
(fr[1]-x[1]*x[2]-x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[1]+
(fr[2]-(1-x[1])*x[2]+x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[2]+
(fr[3]-x[1]*(1-x[2])+x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[3]+
(fr[4]-(1-x[1])*(1-x[2])-x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[4]
}
obj<-function(x){
ifelse(fr[1]==0,0,(fr[1]-x[1]*x[2]-x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[1])+
ifelse(fr[2]==0,0,(fr[2]-(1-x[1])*x[2]+x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[2])+
ifelse(fr[3]==0,0,(fr[3]-x[1]*(1-x[2])+x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[3])+
ifelse(fr[4]==0,0,(fr[4]-(1-x[1])*(1-x[2])-x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[4])
}
sval=rep(0.1,3)
fr<-c(0,0.1,0.2,0.3)
optim(sval,obj, method="BFGS")$par
Thank you,
Serguei
--
___________________________________________________________________
Austrian Institute of Economic Research (WIFO)
Name: Serguei Kaniovski P.O.Box 91
Tel.: +43-1-7982601-231 Arsenal Objekt 20
Fax: +43-1-7989386 1103 Vienna, Austria
Mail: Serguei.Kaniovski at wifo.ac.at
http://www.wifo.ac.at/Serguei.Kaniovski
Have you tried making contour plots of your objective function,
e.g., using expand.grid and 'contour' or 'contourplot', as
described in
Venables and Ripley (2002) Modern Applied Statistics with S (Springer)?
Hope this helps.
Spencer Graves
Serguei Kaniovski wrote:> Hi all,
>
> I apply "optim" to the function "obj", which minimizes
the goodness of
> fit statistic and obtains Pearson minimum chi-squared estimate for x[1],
> x[2] and x[3]. The vector "fr" contains the four observed
frequencies.
>
> Since "fr[i]" appears in the denominator, I would like to
substitute "0"
> in the sum if fr[i]=0.
>
> I tried an "ifelse" condition which works, but in some cases the
> solution seems rather odd. Is there anything wrong with the way second
> objective function is coded?
>
> obj<-function(x){
> (fr[1]-x[1]*x[2]-x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[1]+
> (fr[2]-(1-x[1])*x[2]+x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[2]+
> (fr[3]-x[1]*(1-x[2])+x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[3]+
> (fr[4]-(1-x[1])*(1-x[2])-x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[4]
> }
>
> obj<-function(x){
>
ifelse(fr[1]==0,0,(fr[1]-x[1]*x[2]-x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[1])+
>
ifelse(fr[2]==0,0,(fr[2]-(1-x[1])*x[2]+x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[2])+
>
ifelse(fr[3]==0,0,(fr[3]-x[1]*(1-x[2])+x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[3])+
>
ifelse(fr[4]==0,0,(fr[4]-(1-x[1])*(1-x[2])-x[3]*sqrt(x[1]*(1-x[1])*x[2]*(1-x[2])))^2/fr[4])
> }
>
> sval=rep(0.1,3)
> fr<-c(0,0.1,0.2,0.3)
>
> optim(sval,obj, method="BFGS")$par
>
> Thank you,
> Serguei
>