new; 
cls;

load x[150,45]=c:/data/b1/data1.csv;
output file=c:/data/data1.out reset;

b0=zeros(2,1);

library maxlik;
maxset;
_max_Algorithm=4;
{b,ff,gg,cov,retcode} = maxlik(x,0,&ll,b0);
call maxprt(b,ff,gg,cov,retcode);

/* out put */
output on;

format /rd 8,5;
print "L(0) ="  ll(b0,x);
print "L(b) ="  ll(b,x);
print "ro^2 ="  (ll(b0,x)-ll(b,x))/ll(b0,x);
print "ro-^2="  (ll(b0,x)-(ll(b,x)-rows(b0)))/ll(b0,x);
print;
print  prob(b,x);

output off;

proc prob(b0,x);
    local v,sv;
    local p;
    
    v=zeros(rows(x),15);
    p=zeros(rows(x),15);

    v[., 1]=exp(x[.,16]*b0[1]+x[.,31]*b0[2]).*x[., 1];
    v[., 2]=exp(x[.,17]*b0[1]+x[.,32]*b0[2]).*x[., 2];
    v[., 3]=exp(x[.,18]*b0[1]+x[.,33]*b0[2]).*x[., 3];
    v[., 4]=exp(x[.,19]*b0[1]+x[.,34]*b0[2]).*x[., 4];
    v[., 5]=exp(x[.,20]*b0[1]+x[.,35]*b0[2]).*x[., 5];
    v[., 6]=exp(x[.,21]*b0[1]+x[.,36]*b0[2]).*x[., 6];
    v[., 7]=exp(x[.,22]*b0[1]+x[.,37]*b0[2]).*x[., 7];
    v[., 8]=exp(x[.,23]*b0[1]+x[.,38]*b0[2]).*x[., 8];
    v[., 9]=exp(x[.,24]*b0[1]+x[.,39]*b0[2]).*x[., 9];
    v[.,10]=exp(x[.,25]*b0[1]+x[.,40]*b0[2]).*x[.,10];
    v[.,11]=exp(x[.,26]*b0[1]+x[.,41]*b0[2]).*x[.,11];
    v[.,12]=exp(x[.,27]*b0[1]+x[.,42]*b0[2]).*x[.,12];
    v[.,13]=exp(x[.,28]*b0[1]+x[.,43]*b0[2]).*x[.,13];
    v[.,14]=exp(x[.,29]*b0[1]+x[.,44]*b0[2]).*x[.,14];
    v[.,15]=exp(x[.,30]*b0[1]+x[.,45]*b0[2]).*x[.,15];    
    sv=sumc(v');
    
    p[., 1]=v[., 1]./sv;
    p[., 2]=v[., 2]./sv;
    p[., 3]=v[., 3]./sv;
    p[., 4]=v[., 4]./sv;
    p[., 5]=v[., 5]./sv;
    p[., 6]=v[., 6]./sv;
    p[., 7]=v[., 7]./sv;
    p[., 8]=v[., 8]./sv;
    p[., 9]=v[., 9]./sv;
    p[.,10]=v[.,10]./sv;
    p[.,11]=v[.,11]./sv;
    p[.,12]=v[.,12]./sv;
    p[.,13]=v[.,13]./sv;
    p[.,14]=v[.,14]./sv;
    p[.,15]=v[.,15]./sv;

    retp(p);
endp;

proc ll(b0,x);
    local pp,cc,al;
    cc=ones(rows(x),1)~zeros(rows(x),14);
    pp=prob(b0,x);
    al=sumc(ln(sumc((pp.*x[.,1:15].*cc)')));
    retp(al);
endp;