# Estimating a mixture of probability functions
# for a given set of data (sample) following a probability distribution
# find the parameters to maximize the probabity (likelihood)
#
yed<- read.table("http://web.pdx.edu/~crkl/ec575/data/YED20.TXT",header=T,nrows=20)
# file name may be case sensitive!
# yed<- read.table("C:/course17/ec575/data/yed20.txt",header=T,nrows=20)
summary(yed)
y=yed$Income/10  # scale the data may help, y ~ prob. dist.

# mixture of two normal distributions
# b[1]=mu1, b[2]=sigma1, b[3]=mu2, b[4]=sigma2
llf1<-function(b) {
  p<-pnorm(b[5])
  log(p*dnorm(y,b[1],b[2])+(1-p)*dnorm(y,b[3],b[4]))
}

llf2<-function(b) {
  p<-exp(b[5])/(1+exp(b[5]))   # logist transformation
  log(p*dnorm(y,b[1],b[2])+(1-p)*dnorm(y,b[3],b[4]))
}

library(maxLik)  # need to install maxLik package
M1<-maxLik(llf1,start=c(3,2,2,1,0))
summary(M1)
pnorm(M1$estimate[5])  # lambda

M2<-maxLik(llf2,start=c(3,2,2,1,0))
summary(M2)
exp(M2$estimate[5])/(1+exp(M2$estimate[5]))  # lambda

# optim(c(3,2,2,1,0.5),llf,method="BFGS",hessian=T,control=list(fnscale=-1))