标题: Bayesian Computation with SAS 1. [打印本页] 作者: shiyiming 时间: 2011-1-4 21:05 标题: Bayesian Computation with SAS 1. From oloolo's blog on SasProgramming
<p><a href="http://feedads.g.doubleclick.net/~a/T_SJikcyef9C6jgJlkySU_RAe48/0/da"><img src="http://feedads.g.doubleclick.net/~a/T_SJikcyef9C6jgJlkySU_RAe48/0/di" border="0" ismap="true"></img></a><br/>
<a href="http://feedads.g.doubleclick.net/~a/T_SJikcyef9C6jgJlkySU_RAe48/1/da"><img src="http://feedads.g.doubleclick.net/~a/T_SJikcyef9C6jgJlkySU_RAe48/1/di" border="0" ismap="true"></img></a></p>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 and I am going to do this for next several month to cover all chapters in the book. Readers are encouraged to buy this book to understand what's behind the code. The post here will only spend minimum effects for that purpose.<br />
<br />
This post will cover Chapter 1.<br />
<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;
/* Section 1.2.2 Read Tab Deliminated Data into SAS */
libname data "&datafolder";
proc import datafile="&datafolder\studentdata.txt" out=student dbms=dlm replace;
datarow=2;
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";
/* section 1.2.2 To see Variable Names and print first row of this data */
proc print data=student(obs=1) split="*"; run;
/* section 1.2.3 To summarize and graph a single Batch */
proc freq data=student;
table Drink;
run;
/* section 1.2.3 Barplot on single variable */
proc freq data=student noprint;
table Drink /out=t;
run;
proc gbarline data=t;
bar Drink /discrete space=4 sumvar=count;
run;quit;
/* Actually, in SAS, it is better off to directly use PROC GBARLINE */
proc gbarline data=student;
bar Drink /discrete space=3 ;
run;quit;
/* section 1.2.4 Compare Batch */
data student;
set student;
hours_of_sleep=WakeUp-ToSleep;
label hours_of_sleep="Hours Of Sleep";
run;
proc means data=student min q1 median q3 max nmiss maxdec=2;
var hours_of_sleep;
run;
/* section 1.2.4 Histogram */
proc univariate data=student noprint;
var hours_of_sleep;
histogram /midpoints=2.5 to 12.5 by 1;
run;
proc means data=student
min q1 median q3 max nmiss
nway maxdec=2;
class Gender;
var Haircut;
run;
%macro jitter(var, newname, data=last, factor=1, amount=);
/* follow the JITTER function in R.
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;
/* 1.2.5 Study the relationship between variables */
ods select all ;
proc sgplot data=student;
scatter x=xts y=xhos;
reg x=ToSleep y=Hours_of_sleep;
run;
ods html off;
/* 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;