Bayesian Computation with SAS (1).

shiyiming · 发表于 2011-1-5 00:06:17

From oloolo's blog on SasProgramming

<a href="http://feedads.g.doubleclick.net/~a/Q90LzmIt3V1FHewfGpr5IQhxZV0/0/da"><img src="http://feedads.g.doubleclick.net/~a/Q90LzmIt3V1FHewfGpr5IQhxZV0/0/di" border="0" ismap="true"></img></a> 
<a href="http://feedads.g.doubleclick.net/~a/Q90LzmIt3V1FHewfGpr5IQhxZV0/1/da"><img src="http://feedads.g.doubleclick.net/~a/Q90LzmIt3V1FHewfGpr5IQhxZV0/1/di" border="0" ismap="true"></img></a>The book "<a href="http://bayes.bgsu.edu/bcwr/">Bayesian Computation with R</a>" by <a href="http://bayes.bgsu.edu/">Jim Albert</a> is an easy to read entry level book on applied Bayesian Statistics. While the book was written for R users, it is not difficult to translate the languages between R and SAS and I believe it is a good way to show Bayesian capability of SAS. In the next several months, I am going to translate all R code in the book into SAS. Readers are encouraged to buy this book to understand what's behind the code. The posts here will only spend minimum effects to explain the statistics underlying. 
 
This post will cover Chapter 1. 
 
We first follow section 1.2.2 
<pre style="background-color: #ebebeb; border-bottom: #999999 1px dashed; border-left: #999999 1px dashed; border-right: #999999 1px dashed; border-top: #999999 1px dashed; color: #000001; font-family: Andale Mono, Lucida Console, Monaco, fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding-bottom: 5px; padding-left: 5px; padding-right: 5px; padding-top: 5px; width: 100%;"><code>
/* Chapter 1 of Bayesian Computation with R */

%let folder=C:\Documents\SAS for Bayesian Compu with R;
%let datafolder=&folder\data;

libname data "&datafolder";
options missing='.' formchar = "|----|+|---+=|-/\<>*" errors=2 fullstimer;

/*@@@@@@ Section 1.2.2 @@@@@*/
/*--> Read Tab Deliminated Data into SAS */
proc import datafile="&datafolder/studentdata.txt" out=student dbms=dlm replace;
 datarow=2; /* redundent with GETNAMES= statement */
guessingrows=2; /* Only use first 2 rows of data to guess the type of data */
 delimiter='09'x; /* this is redudent if we specify DBMS=TAB in PROC IMPORT option */
getnames=YES; /* This is the same as 'header=True' in R's table.read function */
run;

ods pdf file="&folder/Chapter1.pdf";
/*--> To see Variable Names and print first row of this data
Use SPLIT="*" to specify the split character being "*", which controls
line breaks in column headings
*/
proc print data=student(obs=1) split="*"; run;
</code></pre> 
SAS output looks like: 
<div class="separator" style="clear: both; text-align: center;"></div> 
 
Section 1.2.3, summarize frequency tables and use Boxplot to visualize the counts: 
<pre style="background-color: #ebebeb; border-bottom: #999999 1px dashed; border-left: #999999 1px dashed; border-right: #999999 1px dashed; border-top: #999999 1px dashed; color: #000001; font-family: Andale Mono, Lucida Console, Monaco, fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding-bottom: 5px; padding-left: 5px; padding-right: 5px; padding-top: 5px; width: 100%;"><code>

/*@@@@ section 1.2.3 @@@*/
/*--> To summarize and graph a single Batch */
proc freq data=student;
 table Drink;
run;

/*--> Barplot on single variable */
title "Original idea, summarize and plot";
proc freq data=student noprint;
 table Drink /out=t;
run;
goptions reset=all noborder;
proc gbarline data=t;
 bar Drink /discrete space=4 sumvar=count;
run;quit;
title;
/*--> Actually, in SAS, it is better off to directly use PROC GBARLINE */
goptions reset=all noborder;
title "In SAS, PROC GBARLINE directly summarize";
proc gbarline data=student;
 bar Drink /discrete space=3 ;
run;quit;
title;

data student;
 set student;
hours_of_sleep=WakeUp-ToSleep;
label hours_of_sleep="Hours Of Sleep"; /* Label can be useful for annotation purpose */
run;

proc means data=student
 min q1 median q3 max nmiss /* Request same statistics as in the book */
 maxdec=2; /* MAXDEC= specifies keeping 2 decimals */
 var hours_of_sleep;
run;

/*--> Histogram using PROC UNIVARIATE */
proc univariate data=student noprint;
 var hours_of_sleep;
histogram /midpoints=2.5 to 12.5 by 1;
run;

proc sgplot data=student;
 histogram hours_of_sleep/scale=count; /* new SGPLOT is a handy way to draw histogram and density */
density hours_of_sleep/scale=count;
run;
</code></pre> 
Compare SAS outputs: 
 
<pre style="background-color: #ebebeb; border-bottom: #999999 1px dashed; border-left: #999999 1px dashed; border-right: #999999 1px dashed; border-top: #999999 1px dashed; color: #000001; font-family: Andale Mono, Lucida Console, Monaco, fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding-bottom: 5px; padding-left: 5px; padding-right: 5px; padding-top: 5px; width: 100%;"><code>

/*@@@ section 1.2.4 @@@@*/
/*---> Compare Batches using Boxplot */
proc sort data=student out=studentsorted; by Gender; run;
proc boxplot data=studentsorted;
 plot hours_of_sleep*Gender;
run;

proc means data=student
 min q1 median q3 max nmiss
 nway maxdec=2;
 class Gender;
 var Haircut;
run;
</code></pre> 
 
 
<pre style="background-color: #ebebeb; border-bottom: #999999 1px dashed; border-left: #999999 1px dashed; border-right: #999999 1px dashed; border-top: #999999 1px dashed; color: #000001; font-family: Andale Mono, Lucida Console, Monaco, fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding-bottom: 5px; padding-left: 5px; padding-right: 5px; padding-top: 5px; width: 100%;"><code>

%macro jitter(var, newname, data=last, factor=1, amount=);
/* follow the JITTER function in R, no fuzzy applied yet .
if amount is given, then use this value for disturbance
else if amount is 0, then use the range of &var for disturbance
else if amount is NULL, then use the smallest difference among
 distinct values of &var for disturbance. if all values are
 the same, then use the value obtained from range

if range of &var is 0, then use the lower range for disturbance
otherwise
*/
%let blank=%str( );
%if &data=' ' %then %let data=last;
%local fid;

data _null_;
 fid=round(ranuni(0)*10000, 1);
 call symput('fid', compress(fid));
run;
proc means data=&data noprint;
 var &var;
 output out=_util&fid range(&var)=range min(&var)=min max(&var)=max;
run;
data _util&fid;
 set _util&fid;
 if range^=0 then z=range; else z=min;
 if z=0 then z=1;
run;

%if %eval(&amount=&blank) %then %do;
proc sort data=&data.(keep=&var where=(&var^=.)) out=_xuti&fid nodupkey;
 by &var;
run;
data _duti&fid;
 set _xuti&fid nobs=ntotal end=eof;
array _x{2} _temporary_;
if ntotal=1 then do;
 amount=&factor/50*abs(&var);
keep amount;
output;
stop;
end;
else do;
 if _n_=1 then do; _x[1]=&var; _x[2]=constant('BIG'); end;
else do;
 _x[2]=min(_x[2], &var - _x[1]);
 _x[1]=&var;
 if eof then do;
 amount=&factor/5*abs(_x[2]);
 keep amount;
 output;
 end;
end;
end;
 run;

%end;
%else %if %eval(&amount=0) %then %do;
data _duti&fid;
 set _util&fid;
amount=&factor*z/50;
keep amount;
output;
run;
%end;
%else %do;
data _duti&fid;
 amount=&amount;
keep amount;
output;
run;
%end;

proc sql noprint;
 select name into :keepvars separated by ' '
 from sashelp.vcolumn
 where libname='WORK' and memname=%upcase("&data") and memtype="DATA"
 ;
quit;
data &data;
 array _x{1} _temporary_;
 if _n_=1 then do;
 set _duti&fid;
 _x[1]=amount;
 end;
 set &data;
 &newname=&var + ranuni(0)*(2*_x[1])-_x[1];
 label &newname="jitter(&var)";
 keep &keepvars &newname;
run;
proc datasets library=work nolist;
 delete _duti&fid _xuti&fid _util&fid;
run;quit;
%mend;
</code></pre> 
 
 
<pre style="background-color: #ebebeb; border-bottom: #999999 1px dashed; border-left: #999999 1px dashed; border-right: #999999 1px dashed; border-top: #999999 1px dashed; color: #000001; font-family: Andale Mono, Lucida Console, Monaco, fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding-bottom: 5px; padding-left: 5px; padding-right: 5px; padding-top: 5px; width: 100%;"><code>

%jitter(hours_of_sleep, xhos, data=student, factor=1, amount=);
%jitter(ToSleep, xts, data=student, factor=1, amount=);

goptions border reset=all;
axis1 label=("jitter(ToSleep)") major=(height=2 number=5) minor=none;
axis2 label=(angle=90 "jitter(Hours of Sleep") major=(height=2 number=5) minor=none;
symbol1 value=circle interpol=none color=black;
proc gplot data=student;
 plot xhos*xts/haxis=axis1 vaxis=axis2;
run;quit;
</code></pre> 
 
 
<pre style="background-color: #ebebeb; border-bottom: #999999 1px dashed; border-left: #999999 1px dashed; border-right: #999999 1px dashed; border-top: #999999 1px dashed; color: #000001; font-family: Andale Mono, Lucida Console, Monaco, fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding-bottom: 5px; padding-left: 5px; padding-right: 5px; padding-top: 5px; width: 100%;"><code>

/* old school */

proc reg data=student lineprinter;
 title "Line Printer Style";
 model Hours_of_Sleep=ToSleep /noprint;
var xhos xts;
plot pred.*xts='X' xhos*xts/overlay symbol='o';
run;quit;
title;

ods select none;
title "ODS Graphics";
ods graphics on;
proc reg data=student ;
 model Hours_of_Sleep=ToSleep;
var xhos xts;
output out=student pred=p;
 ods select FitPlot;
run;quit;
ods graphics off;
title;
ods select all;

proc sort data=student out=studentsorted;
 by ToSleep;
run;
goptions border reset=all;
axis1 label=("jitter(ToSleep)") order=(-2 to 6 by 2) major=(height=2 number=5) minor=none;
axis2 label=(angle=90 "jitter(Hours of Sleep") offset=(, 2 pct) major=(height=2 number=5) minor=none;
symbol1 value=circle interpol=none color=black;
symbol2 value=none interpol=line color=black;
proc gplot data=studentsorted;
 plot xhos*xts p*xts/overlay haxis=axis1 vaxis=axis2;
run;quit;

proc sgplot data=student;
 scatter x=xts y=xhos;
reg x=ToSleep y=Hours_of_Sleep;
run;

</code></pre> 
 
 
<pre style="background-color: #ebebeb; border-bottom: #999999 1px dashed; border-left: #999999 1px dashed; border-right: #999999 1px dashed; border-top: #999999 1px dashed; color: #000001; font-family: Andale Mono, Lucida Console, Monaco, fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding-bottom: 5px; padding-left: 5px; padding-right: 5px; padding-top: 5px; width: 100%;"><code>

/* section 1.3 Robustness of T-Statistics */
/* section 1.3.2 Write a function to compute T-Statistics */
%macro tstatistics(vars, dsn=last, outdsn=, class=);
%if &class='' %then %let classstmt=;
%else %let classstmt=class &class;

%let nvars=%sysfunc(count(&vars, ' '));

%let var1=%scan(&vars, 1, ' ');
%let var2=%scan(&vars, 2, ' ');
proc means data=&dsn noprint nway;
 &classstmt;
 var &vars;
 output out=tdata(where=(_STAT_ in ('STD', 'MEAN', 'N')));
run;
data &outdsn;
 set tdata; by &class _TYPE_;
 retain m n mean1 mean2 std:;
 select(compress(_STAT_));
 when('N') do;
 m=&var1;
 n=&var2;
end;
when('MEAN') do;
 mean1=&var1;
 mean2=&var2;
end;
 when('STD') do;
 std1=&var1;
 std2=&var2;
end;
otherwise;
 end;
 if last._TYPE_ then do;
 sp=sqrt(((m-1)*std1**2 + (n-1)*std2**2)/(m+n-2));
 t= (mean1-mean2)/(sp * sqrt(1/m + 1/n));
 keep &class _TYPE_ m n std: mean: sp t;
 output;
 end;
run;
%mend;

filename BCR "&folder";
%inc BCR(tstatistics);

data test;
 input x y;
cards;
1 5
4 4
3 7
6 6
5 10
;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=);

/* section 1.3.3 Monte Carlo study on the Robustness of T-Statistics */
%let alpha=0.1;
%let m=10;
%let n=10;
%let Niter=10000;
data test;
 retain seed1 99787 seed2 76568;
do iter=1 to &Niter;
 do j=1 to max(&n, &m);
 if j<=&m then call rannor(seed1, x); else x=.;
if j<=&n then call rannor(seed2, y); else y=.;
keep iter x y;
output;
end;
end;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=iter);

data _null_;
 set tdata2 end=eof;
retain n_reject 0;
if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
if eof then do;
 call symput('n_reject', n_reject/_n_);
end;
run;
%put &n_reject;

/* section 1.3.4 Study the robustness of T-Statistics */

/* 1.) both are standard normal */

%let alpha=0.1;
%let m=10;
%let n=10;
%let Niter=10000;
data test;
 retain seed1 987687 seed2 76568;
do iter=1 to &Niter;
 do j=1 to max(&n, &m);
 if j<=&m then call rannor(seed1, x); else x=.;
if j<=&n then call rannor(seed2, y); else y=.;
keep iter x y;
output;
end;
end;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=iter);

data _null_;
 set tdata2 end=eof;
retain n_reject 0;
if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
if eof then do;
 call symput('n_reject', n_reject/_n_);
end;
run;
%put &n_reject;

/* 2.) both are normal with std=1 and 10 */

%let alpha=0.1;
%let m=10;
%let n=10;
%let Niter=10000;
data test;
 retain seed1 99787 seed2 76568;
do iter=1 to &Niter;
 do j=1 to max(&n, &m);
 if j<=&m then call rannor(seed1, x); else x=.;
if j<=&n then call rannor(seed2, y); else y=.;
y=y*10;
keep iter x y;
output;
end;
end;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=iter);

data _null_;
 set tdata2 end=eof;
retain n_reject 0;
if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
if eof then do;
 call symput('n_reject', n_reject/_n_);
end;
run;
%put &n_reject;

/* 2.) both are normal with std=1 and 10 */

%let alpha=0.1;
%let m=10;
%let n=10;
%let Niter=10000;
data test;
 retain seed1 99787 seed2 76568;
do iter=1 to &Niter;
 do j=1 to max(&n, &m);
 if j<=&m then call rannor(seed1, x); else x=.;
if j<=&n then call rannor(seed2, y); else y=.;
y=y*10;
keep iter x y;
output;
end;
end;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=iter);

data _null_;
 set tdata2 end=eof;
retain n_reject 0;
if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
if eof then do;
 call symput('n_reject', n_reject/_n_);
end;
run;
%put &n_reject;

/* 3.) both are T-distribution with df=4 */

%let alpha=0.1;
%let m=10;
%let n=10;
%let Niter=10000;
data test;
 retain seed1 99787 seed2 76568;
call streaminit(seed1);
do iter=1 to &Niter;
 do j=1 to max(&n, &m);
 if j<=&m then x=rand('T', 4); else x=.;
if j<=&n then y=rand('T', 4); else y=.;
keep iter x y;
output;
end;
end;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=iter);

data _null_;
 set tdata2 end=eof;
retain n_reject 0;
if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
if eof then do;
 call symput('n_reject', n_reject/_n_);
end;
run;
%put &n_reject;

/* 4.) one normal with mu=10 and std=2 while the other is exponential with rate=0.10 */

%let alpha=0.1;
%let m=10;
%let n=10;
%let Niter=10000;
data test;
 retain seed1 99787 seed2 76568;
do iter=1 to &Niter;
 do j=1 to max(&n, &m);
 if j<=&m then call rannor(seed1, x); else x=.;
if j<=&n then call ranexp(seed2, y); else y=.;
x=10+x*2;
y=y/0.1;
keep iter x y;
output;
end;
end;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=iter);

data _null_;
 set tdata2 end=eof;
retain n_reject 0;
if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
if eof then do;
 call symput('n_reject', n_reject/_n_);
end;
run;
%put &n_reject;

/* draw the overlay density curves of empirical T-statistics and T-distribution */
proc sort data=tdata2;
 by t;
run;
data tdata2;
 set tdata2;
dt=pdf('T', t, 18);
run;

ods select none;
ods output Controls=Controls;
ods output UnivariateStatistics=UniStat;
proc kde data=tdata2 ;
 univar t /out=density method=snr unistats;
run;
ods select all;

data density;
 set density;
dt=pdf('T', value, 18);
run;

data _null_;
 set UniStat;
if Descr='Bandwidth' then call symput ('bw', compress(round(t, 0.0001)));
run;

goptions reset=all border;

title "Comparing Empirical and Theoretical Densities";
legend label=none position=(top right inside) mode=share;
axis1 label=("N=&Niter Bandwidth=&bw") order=(-4 to 8 by 2) offset=(1,1) major=(height=2) minor=none;
axis2 label=(angle=90 "Density") order=(0 to 0.4 by 0.1) offset=(0, 0.1) major=(height=2) minor=none;
symbol1 interpol=join color=black value=none width=3;
symbol2 interpol=join color=grey value=none width=1;
proc gplot data=density;
 plot density*value dt*value/overlay haxis=axis1 vaxis=axis2 legend=legend;
 label dt='t(18)' density='Exact';
run;quit;

ods pdf close;
</code></pre> 
<a href="http://www.amazon.com/Bayesian-Computation-R-Use/dp/0387922970?ie=UTF8&tag=xie1978&link_code=bil&camp=213689&creative=392969" imageanchor="1" target="_blank"><img alt="Bayesian Computation with R (Use R)" src="http://ws.amazon.com/widgets/q?MarketPlace=US&ServiceVersion=20070822&ID=AsinImage&WS=1&Format=_SL160_&ASIN=0387922970&tag=xie1978" /></a><img alt="" border="0" height="1" src="http://www.assoc-amazon.com/e/ir?t=xie1978&l=bil&camp=213689&creative=392969&o=1&a=0387922970" style="border-bottom: medium none; border-left: medium none; border-right: medium none; border-top: medium none; margin: 0px; padding-bottom: 0px! important; padding-left: 0px! important; padding-right: 0px! important; padding-top: 0px! important;" width="1" /><div class="blogger-post-footer"><img width='1' height='1' src='https://blogger.googleusercontent.com/tracker/29815492-6025913615236309735?l=www.sas-programming.com' alt='' /></div><img src="http://feeds.feedburner.com/~r/SasProgramming/~4/yL6bIedyqD4" height="1" width="1"/>

		自动登录	找回密码
密码			立即注册