data biossay;
input doses rats deaths freq;
y=log(freq/(1-freq)); x=doses;
cards;
-0.863    5    0.5  0.1
-0.296    5    1.0  0.2
-0.053    5    3.0  0.6
0.727     5    4.5  0.9
run;

ods output ParameterEstimates = b;
proc reg data=biossay;
model y=x;
run;quit;


proc iml;
use biossay; 
read all into biossay;
use b;
read all var{Estimate,StdErr} into parameter;

start input (LL) global(alpha, beta);
alpha = do(-2,2,(4/(LL-1)) ) ;
beta  = do(-5,10,(15/(LL-1)));
finish input;

start log_post(alpha,beta,dd);
  ldens = 0; 
  do i=1 to nrow(dd);
    theta = exp(alpha+beta*dd[i,1])/(1+exp(alpha+beta*dd[i,1]));
    ldens = ldens + dd[i,3]*log(theta)+(dd[i,2]-dd[i,3])*log(1-theta);
   end;
  return (ldens);
finish; 

start grid_value(DD,alpha,beta);
  ll = ncol(alpha);
  PP =J(ll,ll,.);
  do i=1 to ll; 
    do j = 1 to ll;
      PP[i,j]= exp(log_post(alpha[i],beta[j],dd));	
    end;
 end;
 return(PP/sum(PP));
finish;

start sampling(M,PP,alpha,beta) global(post_alpha,post_beta);
  alpha_mar = PP[,+];
  alpha_cdf = cusum(alpha_mar); 
  post_alpha = J(M,1,0);
  post_beta  = J(M,1,0);
  do i = 1 to M;
    uuu = ranuni(0);
    Fhat_alpha  =  max( alpha_cdf #(alpha_cdf <= uuu)); 
    post_alpha[i] = alpha[loc (alpha_cdf = Fhat_alpha)]; 
 
    junk  = nrow(alpha[loc (alpha <= post_alpha[i,1]) ]);
    PP[junk, ]  =  PP[junk,]/PP[junk,+];                                   
    beta_cond_cdf =   cusum(PP[junk,]);
    uuu = ranuni(0);
    Fhat_beta  =  max( beta_cond_cdf[loc(beta_cond_cdf < uuu)]);    
    post_beta[i] = beta[loc (beta_cond_cdf = Fhat_beta)]; 
  end;
finish; 

start joint_post(alpha,beta,DD);
k1= ncol(alpha);k2=ncol(beta); logdens=j(k1*k2,1,0);
k=0; 
do i=1 to k1;
   do j=  1 to k2;
    k=k+1;   logdens[k]  = log_post(alpha[i],beta[j],DD); 
   end;
end;
dens = exp(logdens-max(logdens)); 
contour_data=t(alpha)@J(k2,1,1)||J(k1,1,1)@t(beta)||dens;
return (contour_data);
finish;  

run input(200);
PP=grid_value(biossay, alpha, beta);
run  sampling(1000,PP,alpha,beta);
X= joint_post(alpha,beta, biossay);

create work.contour var{alpha, beta, density};
append from X; free X;

X=post_alpha||post_beta;

create work.posterior var{alpha, beta};
append from X; free /; 
quit;

data annopoints;                                                                                                                        
  set posterior;                                                                                                                             
  function='symbol';
  text='circle';
  xsys='2'; ysys='2'; zsys='2';
  x=alpha; y=beta;
  when='a';  size=0.2;                                                                                                             
run;         
proc gcontour data=contour;
  title " Posterior Density ";
   axis1 label=("Alpha ") order=(-2 to 2 by 1) offset=(0) minor=none ;
   axis2 label=("Beta ") order=(-5 to 10 by 5)  minor=none;  
   plot beta*alpha=density/nlevels=10 haxis=axis1 vaxis=axis2 noframe autolabel=(check=none) nolegend  caxis=black annotate=annopoints;
run; title " "; quit;