/*
** Lesson 6.4: Mixture of Two Normal Distributions
** See Greene [1999], Chapter 4
*/
use gpe2;
output file=output6.4 reset;
load data[21,2]=gpe\yed20.txt;
x=data[2:21,1]/10; @ income data: scaling may help @

call reset;
_nlopt=1;
_method=5;
_iter=100;
_b={3,3,2,2,0.5};
call estimate(&llf,x);
p1=prob(x,__b);

_b={3,3,2,2,0};
call estimate(&llfn,x);
p2=probn(x,__b);

print p1~p2;

end;

/*
mixture of two normal distributions
mu1=b[1], mu2=b[2]
se1=b[3], se2=b[4]
prob.(drawn from the 1st distribution)=b[5]
*/
proc llf(x,b);
    local pdf1,pdf2;
    pdf1=pdfn((x-b[1])/b[3])/b[3];
    pdf2=pdfn((x-b[2])/b[4])/b[4];
    retp(sumc(ln(b[5]*pdf1+(1-b[5])*pdf2)));
endp;

/*
using cdfn transformation with unrestricted b[5]
prob(R=1) = prob(R* = b[5] + e >0), and
prob(R=0) = prob(R* = b[5] + e <=0 ), e is n(0,1)
prob(e|e>-b[5]) = cdfn(b[5]) = prob(R=1)
prob(e|e<=-b[5]) = 1-cdfn(b[5]) = prob(R=0)
*/
proc llfn(x,b);
    local pdf1,pdf2;
    pdf1=pdfn((x-b[1])/b[3])/b[3];
    pdf2=pdfn((x-b[2])/b[4])/b[4];
    retp(sumc(ln(cdfn(b[5])*pdf1+(1-cdfn(b[5]))*pdf2)));
endp;

proc prob(x,b);
    local pdf1,pdf2,pdf;
    pdf1=pdfn((x-b[1])/b[3])/b[3];
    pdf2=pdfn((x-b[2])/b[4])/b[4];
    pdf=b[5]*pdf1+(1-b[5])*pdf2;
    retp((b[5]*pdf1./pdf)~((1-b[5])*pdf2./pdf));
endp;

proc probn(x,b);
    local pdf1,pdf2,pdf;
    pdf1=pdfn((x-b[1])/b[3])/b[3];
    pdf2=pdfn((x-b[2])/b[4])/b[4];
    pdf=cdfn(b[5])*pdf1+(1-cdfn(b[5]))*pdf2;
    retp((cdfn(b[5])*pdf1./pdf)~((1-cdfn(b[5]))*pdf2./pdf));
endp;