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

/* normal probability distribution: b[1]=mu, b[2]=sigma */
fn llfn(x,b)=sumc(ln(pdfn((x-b[1])/b[2])./b[2]));
fn llfln(x,b)=sumc(ln(pdfn((ln(x)-b[1])/b[2])./(b[2].*x)));

/* gamma probability distribution: b[1]=rho, b[2]=lambda */
fn pdfg(x,b)=((b[2]^b[1])./gamma(b[1])).*exp(-b[2]*x).*x^(b[1]-1);
fn llfg(x,b)=sumc(ln(pdfg(x,b)));

call reset;
_nlopt=1;
_method=4;
_iter=100;

_b={3.0,2.0};
call estimate(&llfn,x);

_b={1.0,0.5};
call estimate(&llfln,x);

_b={2.0,0.5};
call estimate(&llfg,x);
end;