Bayesian Computation with SAS (2).

shiyiming

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&gt;0 &amp; p&lt;1 then do;
      log_like=log(p)*&amp;n1 + log(1-p)*&amp;n0;   
   end;
   else
            log_like=-999*(p=0)*(&amp;n1&gt;0) + (p=1)*(&amp;n0&gt;0);   
   log_posterior=log(prior) + log_like;
   if log_posterior&gt;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;

symbol value=NONE interpol=Needle;
axis1  label=(angle=90  'Posterior Probability');
axis2  label=('p');
proc gplot data=p2x;
      plot  posterior*p /vaxis=axis1  haxis=axis2;
run;quit;


/* beta prior */

%let a=3.4;
%let b=7.4;
%let s=11;
%let f=16;

data p;
      do i=1 to 500;
      p=i/500;
   prior=pdf('beta', p, &amp;a, &amp;b);
   like=pdf('beta', p, &amp;s, &amp;f);
   post=pdf('beta', p, &amp;a+&amp;s, &amp;b+&amp;f);
   drop i;
   output;
   end;
run;

goptions reset=all;
symbol1 interpol=j  color=red  width=1  value=none;
symbol2 interpol=j  color=blue width=2  value=none;
symbol3 interpol=j  color=gree width=3  value=none;
axis1   label=(angle=90  'Density')  order=(0 to 5 by 1) minor=none;
axis2   label=('p')      order=(0 to 1 by 0.2)  minor=none;;
legend label=none   position=(top right inside)  mode=share;
proc gplot data=p;
      plot post*p prior*p  like*p  /overlay  vaxis=axis1  haxis=axis2  legend=legend;
run;quit;

/* generate 1000 random obs from posterior distribution */
data rs;
      call streaminit(123456);
      do i=1 to 1000;
         x=rand('beta', &amp;a+&amp;s, &amp;b+&amp;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;

proc sgplot data=rs;
      histogram x;
   title "Histogram of ps";
run;
title;
</code></pre><div class="blogger-post-footer"><img width='1' height='1' src='https://blogger.googleusercontent.com/tracker/29815492-4705248398300447375?l=www.sas-programming.com' alt='' /></div><img src="http://feeds.feedburner.com/~r/SasProgramming/~4/1wEPb0UVKbE" height="1" width="1"/>