[R] pb with optimimization and fitdistr

[Ravi Varadhan]

I think the problem lies in the mixture model.  First, why do you have two
free parameters p and q.  Shouldn’t that be p and (1-p)?  Secondly, I am not
sure that your data, d2, is compatible with a binary mixture model. It seems
like a sensible binary mixture model cannot be fitted for your data. It may
be over parametrized.  A model with just a single beta distribution seems to
do ok.  


loglik2 <- function(par, data) {
sh11 <- par[1]
sh12 <- par[2]
sum(dbeta(data,shape1=sh11,shape2=sh12, log=TRUE))
}

optim(par=c(0.1,0.1), fn=loglik2, data=d2, lower=c(0,0),
control=list(fnscale=-1))

Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml



----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Nicolas Turenne
Sent: Wednesday, November 04, 2009 12:17 PM
To: r-help at r-project.org
Subject: [R] pb with optimimization and fitdistr

Hello

i try to fit a data series (N below) with a model consisting of a 
mixture of two beta distributions
for that i am using fitdistr of package MASS

as follows

 > library(MASS)
 > N=c(796,3586,4089,3364,2745,1992,1120,432,99,10,0,0)
 >  d2 = (N-min(N)+0.000001)/(max(N)-min(N)+0.002)
 >  mixtBeta <- function(x,p,q,sh11,sh12,sh21,sh22) 
p*dbeta(x,shape1=sh11,shape2=sh12)+(q)*dbeta(x,shape1=sh21,shape2=sh22)
 >  ff= fitdistr(d2,mixtBeta , 
start=list(p=0.2,q=0.1,sh11=2,sh12=9,sh21=6,sh22=10))
Il y a eu 50 avis ou plus (utilisez warnings() pour voir les 50 premiers)
 > ff
        p              q             sh11           sh12           
sh21           sh22
  2.265173e+03   2.751748e+03   6.823009e-02   1.230362e-01   
7.130326e-01    1.069649e+00

but manually i find another set of parameters
p=0.26;q=0.12; sh11=2.1, sh12=9, sh21=5.8, sh22=9

S=c(0,1,2,3,4,5,6,7,8,9,10,11)/12
 N=c(796,3586,4089,3364,2745,1992,1120,432,99,10,0,0)
 d2 = (N-min(N)+0.001)/(max(N)-min(N)+0.002)
 x = c(1:100)/100
 p=0.26;q=0.12; y = p*dbeta( x, 2.1, 9  ) + q*dbeta( x, 5.8, 9  );
 plot(S,d2,xlim=c(0,1), ylim=c(0,1))
 lines(x, y )

Is there a problem with fitdistr if one uses a function with parameters ?
or a problem of initialization ?

thanks for help to understand.

regards
Nicolas


-- 
Nicolas Turenne - INRA, Unité Mathématique Informatique et Génome UR1077,
F-78350 Jouy-en-Josas
http://genome.jouy.inra.fr/~turenne/nt.html

______________________________________________
R-help at r-project.org 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.