Re: [R] Fit a mixture of lognormal and normal distributions

[Marc Girondot]

Le 18/03/13 16:24, To . . a écrit :

x.1<-rnorm(6000, 2.4, 0.6)x.2<-rlnorm(10000, 1.3,0.1)X<-c(x.1, x.2)
hist(X,100,freq=FALSE, ylim=c(0,1.5))lines(density(x.1), lty=2, 
lwd=2)lines(density(x.2), lty=2, lwd=2)lines(density(X), lty=4)

Here is a solution:

x.1<-rnorm(6000, 2.4, 0.6)
x.2<-rlnorm(10000, 1.3,0.1)

X<-c(x.1, x.2)
hist(X,100,freq=FALSE, ylim=c(0,1.5))
lines(density(x.1), lty=2, lwd=2)
lines(density(x.2), lty=2, lwd=2)
lines(density(X), lty=4)

fitnormlnorm<-function(par, val) {
  p <- 1/(1+exp(-par[5]))
  return(-sum(log(p*dnorm(val, par[1], abs(par[2]), log = FALSE)+
                    (1-p)*dlnorm(val, par[3], abs(par[4]), log = FALSE))))
}

# Mean 1
m1=2.3; s1=0.5
# Mean 2
m2=1.3; s2=0.1
# proportion of 1 - logit transform
p=0

par<-c(m1, s1, m2, s2, p)

result2<-optim(par, fitnormlnorm, method="BFGS", val=X,
               hessian=FALSE, control=list(trace=1))

lines(seq(from=0, to=5, length=100),
      dnorm(seq(from=0, to=5, length=100),
            result2$par[1], result2$par[2]), col="red")

lines(seq(from=0, to=5, length=100),
      dlnorm(seq(from=0, to=5, length=100),
             result2$par[3], result2$par[4]), col="green")

p <- 1/(1+exp(-result2$par[5]))

paste("Proportion of Gaussian data",  p)

lines(seq(from=0, to=5, length=100),
      p*dnorm(seq(from=0, to=5, length=100),
             result2$par[1], result2$par[2])+
        (1-p)*dlnorm(seq(from=0, to=5, length=100),
                    result2$par[3], result2$par[4]), col="blue")

Sincerely

Marc Girondot

--
__________________________________________________________
Marc Girondot, Pr

Laboratoire Ecologie, Systématique et Evolution
Equipe de Conservation des Populations et des Communautés
CNRS, AgroParisTech et Université Paris-Sud 11 , UMR 8079
Bâtiment 362
91405 Orsay Cedex, France

Tel:  33 1 (0)1.69.15.72.30   Fax: 33 1 (0)1.69.15.73.53
e-mail: marc.giron...@u-psud.fr
Web: http://www.ese.u-psud.fr/epc/conservation/Marc.html
Skype: girondot

______________________________________________
R-help@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.