new;

load x[287,23]= c:\sample\kaimono\kaimono.csv;
     x= delif(x,x[.,9].==5208); @pal_k@
     x= delif(x,x[.,9].==5206); @pal_h@
     x= delif(x,x[.,9].==5204); @fuji_s@


library maxlik;
#include maxlik.ext;
maxset;
    

     output file = c:\sample\kaimono\mnl.out reset;
   __output = 2;

clearg prob,selut,sy,ut,jyasuko,takashi,mituko,f_matu,joepla,over_40,woman,sunday,walk,neet,ptj,kaisya,car;

  over_40 = ((x[.,3]).>=40);
  woman   = ((x[.,2]).==2);
  sunday  = ((x[.,6]).==1);
  walk    = ((x[.,23]).==420);
  tuesday = ((x[.,6]).==3);
  neet    = ((x[.,4]).==6);
  ptj     = ((x[.,4]).==5);
  kaisya  = ((x[.,4]).==1);
  car     = ((x[.,23]).==100);
  syufu   = ((x[.,4]).==3);

  
@jyasuko=5272@
@takashi=5212@
@mituko=5213@
@f_matu=5202@
@joepla=5303@

proc  li(b,x);

              @const@    @distance@          @sex@           @age@           
  takashi  = exp(  	     b[5]*x[.,15]/1000                  	         + b[8]*sunday   +b[9]*walk  +b[13]*kaisya );
  mituko   = exp( b[1] + b[5]*x[.,16]/1000                + b[7]*over_40 + b[9]*walk     +b[11]*neet @+b[14]*car@ ) ;  
  f_matu   = exp( b[2] + b[5]*x[.,17]/1000                               + b[8]*sunday             );        
  jyasuko  = exp( b[3] + b[5]*x[.,21]/1000                               + b[10]*tuesday +b[12]*ptj  ); 
  joepla   = exp( b[4] + b[5]*x[.,22]/1000 + b[6]*woman   + b[7]*over_40 );     

  selut= (x[.,9].==5272).*jyasuko + (x[.,9].==5212).*takashi + (x[.,9].==5213).*mituko + (x[.,9].==5202).*f_matu + (x[.,9].==5303).*joepla;
  ut=jyasuko+takashi+mituko+f_matu+joepla;

prob = selut./ut;

retp(ln(prob));

endp;

  start = zeros(13,1);

 let _max_parnames =
@1@ "1__const"
@2@ "2__const"
@3@ "3__const"
@4@ "4__const"
@5@ "distant "
@6@ "sex__wo "
@7@ "age__40 "
@8@ "sunday  "
@9@ "walk  "
@10@"tuesday"
@11@"neet "
@12@"ptj  "
@13@"kaisya "
@@14@"car  "@;

   {b,ff,gg,cov,retcode} = maxlik(x,0,&li,start);
          call maxprt(b,ff,gg,cov,retcode);

_L0 = li(start,x);
L0 = sumc(_L0);

_Lb = li(b,x);
Lb = sumc(_Lb);

print "L(0) =   ";;L0;
print "L(b) =   ";;Lb;
print "ƒÏ^2 =   ";;(L0-Lb)/L0;
print "ƒÏ-^2=   ";;(L0-(Lb-rows(start)))/L0;

@kazoeage@
print "p_takashi=";;sumc(takashi./ut)/266;
print "p_mituko=";;sumc(mituko./ut)/266;
print "p_f_matu=";;sumc(f_matu./ut)/266;
print "p_jyasuko=";;sumc(jyasuko./ut)/266;
print "p_joepla=";;sumc(joepla./ut)/266;

output off;
end;