20150103ブログ用1

conditional probit.sas

proc import out = modelsample
                        datafile = "C:\Users\User\Desktop\三田論\mitaron_conditional_probit1113.csv"
                        dbms=csv REPLACE;
		  getnames=yes;
        datarow=2; 
		run;
proc print data=modelsample(obs=30);
run:
data model(keep=id week year month day brand uicp sicp qicp kicp picp  price cm sns price_rihanari cm_rihanari sns_rihanari choice chosen demo_icp  age sex marry  income1 income2 );
     /*切片*/
     array uicps[5] uicp1 uicp2 uicp3 uicp4 uicp5;
	 array sicps[5] sicp1 sicp2 sicp3 sicp4 sicp5;
	 array qicps[5] qicp1 qicp2 qicp3 qicp4 qicp5;
	 array kicps[5] kicp1 kicp2 kicp3 kicp4 kicp5;
	 array picps[5] picp1 picp2 picp3 picp4 picp5;
	 /*離反なし*/
     array prices[5] uprice sprice qprice kprice pprice;
	 array cms[5] ucm scm qcm kcm pcm;
	 array snss[5] usns ssns qsns ksns psns;
	 /*離反あり*/
     array prices_rihanari[5] uprice_rihanari sprice_rihanari qprice_rihanari kprice_rihanari pprice_rihanari;
	 array cms_rihanari[5] ucm_rihanari scm_rihanari qcm_rihanari kcm_rihanari pcm_rihanari;
	 array snss_rihanari[5] usns_rihanari ssns_rihanari qsns_rihanari ksns_rihanari psns_rihanari;
	 /**/
	 array allbrands[5] $ _temporary_ ('uni' 'shiseido' 'qu' 'kaou' 'pg');
	 set modelsample;
	 do i = 1 to 5;
	     brand = allbrands[i];
		 uicp = uicps[i];
		 sicp = sicps[i];
		 qicp = qicps[i];
		 kicp = kicps[i];
		 picp = picps[i];
		 price = prices[i];
		 cm = cms[i];
		 sns = snss[i];
		 price_rihanari = prices_rihanari[i];
		 cm_rihanari = cms_rihanari[i];
		 sns_rihanari = snss_rihanari[i];
		 choice = (chosen eq brand);
		 output;
	end;
run;
proc print data=model(obs=30);
run;
proc export 
  data=model
  dbms=xlsx 
  outfile="mprobit_ex.xlsx" 
  replace;
run;
/*ods html body='ttest.htm' style=HTMLBlue file="C:\Users\User\Desktop\三田論\mitaron_conditional_probit_ex.html";*/
ods graphics on;
ods pdf  file="C:\Users\Us0er\Desktop\三田論\mitaron_conditional_probit1114.pdf";
proc bchoice
data=model
outpost=mitaron
seed=71
nthreads=8
nbi=30000
nmc=20000
thin=1
hitprob
diag= all
PLOTS=all
stats=all;
class id week year month day chosen(ref='uni') / param=ref order=data; 
model choice = sicp qicp kicp picp /choiceset=(id week year month day chosen ) covprior=iwishart(df=20 scale=100) type=probit;
/*preddist outpred=Predout;*/
random /*sicp qicp kicp picp*/ price cm sns price_rihanari cm_rihanari sns_rihanari / sub=id remean=(/*demo_icp*/ age sex marry income1 income2 /*switch*/) covprior=iwishart(df=20 scale=100) type=un;  
run;
/*%postsum(data=predout)*/
ods pdf close;
/*ods html close;*/
ods graphics off;
proc export
data=mitaron
dbms=csv 
  outfile="mitaron_postsample_ex.csv" 
  replace;
run;