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\mnl3.out reset;
   __output = 2;

clearg prob,selut,sy,ut,jyasuko,takashi,mituko,f_matu,joepla,over_40,woman,tue,car,park,park_car,park2,pub_tran;

  over_40 = x[.,3].>=40;
  woman   = x[.,2].==2;
  tue     = x[.,6].==3;
  park    = x[.,14].==1 .and x[.,23].==420;
  park_car= x[.,14].==2 .and x[.,23].==100;
  park2   = x[.,14].==1;
  pub_tran= x[.,23].==200 .or x[.,23].==240;

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

proc  li(b,x);

              @const@           @distance@          @sex@           @age@        @tue@        @Park@ 
  takashi = exp( b[1] + b[5]*x[.,15]/1000                                         + b[10]                   );  
  mituko  = exp( b[2] + b[5]*x[.,16]/1000               + b[8]*over_40            + b[10]                        );
  f_matu  = exp( b[3] + b[6]*x[.,17]/1000                                                        @+ b[15]*park_car@  );                
  jyasuko = exp( b[4] + b[6]*x[.,21]/1000                                +b[9]*tue                               ); 
  joepla  = exp(        b[6]*x[.,22]/1000 + b[7]*woman  + b[8]*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(10,1);

 let _max_parnames =
@1@ "1__const"
@2@ "2__const"
@3@ "3__const"
@4@ "4__const"
@5@ "dis___1"
@6@ "dis___2"

/*
@5@ "dis___ta"
@6@ "dis___mi"
@7@ "dis___f_"
@8@ "dis___jy"
@9@ "dis___jo"
*/
@10@"sex__wo "
@11@"age_40_t"
@12@"age_40_j"
@13@"tue___jy"
@14@"park__ta"
@15"park__f_"
16"5__const"@
;

   {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;


output off;
end;