标题: Bayesian Computation with SAS (2). [打印本页] 作者: shiyiming 时间: 2011-2-1 23:56 标题: Bayesian Computation with SAS (2). From oloolo's blog on SasProgramming
<p><a href="http://feedads.g.doubleclick.net/~a/PfDDPpDbsupbS4F9hj-QlII4S0g/0/da"><img src="http://feedads.g.doubleclick.net/~a/PfDDPpDbsupbS4F9hj-QlII4S0g/0/di" border="0" ismap="true"></img></a><br/>
<a href="http://feedads.g.doubleclick.net/~a/PfDDPpDbsupbS4F9hj-QlII4S0g/1/da"><img src="http://feedads.g.doubleclick.net/~a/PfDDPpDbsupbS4F9hj-QlII4S0g/1/di" border="0" ismap="true"></img></a></p>In Chapter of the book "Bayesian Computation with R" by Jim Albert, the philosoph of Bayesian statistics is introduced using an example where a parameter regarding proportion of college students that have more than 8 hours sleep a day is studied.<br />
<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 2 of Bayesian Computation with SAS */
/* In this chapter, Jim talked about applying
Bayesian Rule to learn posterior probability
using an example where the proportion of
students sleeping over 8 hours a day is studies.
Bayes Rule:
Posterior(p|data) \prop Prior(p)* likelihood(data|p)
*/
data prior;
input prior @@;
cards;
2 4 8 8 4 2 1 1 1 1
;
run;
data p;
do p=0.05 to 0.95 by 0.1;
output;
end;
run;
data p;
retain sum_prior 0;
do until (eof1);
set prior end=eof1;
sum_prior+prior;
end;
do until (eof2);
merge p prior end=eof2;
prior=prior/sum_prior;
drop sum_prior;
output;
end;
run;
data px;
set p;
output;
prior=0;
output;
run;
goptions reset=all;
symbol interpol=Needle value=none line=1 color=black;
axis1 label=(angle=90 'Prior Probability') order=(0 to 0.3 by 0.1) major=(height=2) minor=none;
axis2 label=('p') order=(0 to 1 by 0.2) minor=none;
proc gplot data=px;
plot prior*p /vaxis=axis1 haxis=axis2;
run;quit;
%let n0=16;
%let n1=11;
data p;
set p end=eof1 nobs=ntotal;
if p>0 & p<1 then do;
log_like=log(p)*&n1 + log(1-p)*&n0;
end;
else
log_like=-999*(p=0)*(&n1>0) + (p=1)*(&n0>0);
log_posterior=log(prior) + log_like;
if log_posterior>max_like then max_like=log_posterior;
run;
proc means data=p noprint;
var log_posterior;
output out=_max max(log_posterior)=max;
run;
data p2;
set _max;
sum_post=0;
do until (eof);
set p end=eof;
sum_post+exp(log_posterior-max);
end;
do until (eof2);
set p end=eof2;
posterior=exp(log_posterior-max)/sum_post;
output;
keep p prior posterior;
end;
run;
data p2x;
set p2;
output;
posterior=0;
output;
run;
/* generate 1000 random obs from posterior distribution */
data rs;
call streaminit(123456);
do i=1 to 1000;
x=rand('beta', &a+&s, &b+&f);
output;
end;
run;
proc univariate data=rs noprint;
var x;
histogram /midpoints=0 to 0.8 by 0.05 cbarline=blue cfill=white outhistogram=hist;
title "Histogram of ps";
run;
title;
axis1 label=(angle=90 "Frequency");
proc gbarline data=hist;
bar _midpt_/discrete sumvar=_count_ space=0 axis=axis1 ;
title "Histogram of ps";
run;
title;