/*                                               */
 /* chex131.sas                                   */
 /* Bayesian Analysis of variance components      */
 /* using Gibbs sampling -- see Example 10.3      */
 /*                                               */
  options ls = 80 ;           * easier to read output ;
  data a ;
   keep f g b t1-t5 t ;
   label f = 'error var'
         g = 'tmt var'
         b = 'mean' ;
   seed = 5151917 ;
 /* set up arrays for vectors n, theta, yb        */
   array yb{5} yb1-yb5 ;
   array  nn{5} nn1-nn5 ;
   array theta{5} t1-t5 ;
 /* here's the data                               */
   p = 5 ;         * number of batches ;
 * b =  ;          * between sum of squares ;
   w = 13.10 ;     * within sum of squares ;
 /* read the rest of the data                     */
   input yb1-yb5 nn1-nn5 ;
 /* prior parameters                              */
   a1 = 1 ;  b1 = 1 ;         * for f = error variance ;
   a2 = 4 ;  b2 = 4 ;         * for g = treatment variance ;
   beta0 = 8 ; phi0 = 16 ;    * for beta = mean ;
 /* starting values for Gibbs sampling            */
   b = 6    ;  f = 5   ; g = 1 ;
   do i = 1 to p ;
   theta(i) = yb(i) ;
   end ;
   put _all_ ;
 /*                                               */
 /* observation loop over parameter values        */
 /*                                               */
   do t = 1 to 2500;     * loop on t ;
 /*                                               */
 /* new parameters of dist of phi given others    */
   a1new = a1 + sum( of nn{*} )/2 ;
   b1new = b1 + w/2 ;
   do i = 1 to p ;
   b1new = b1new + nn(i)*(yb(i)-theta(i))*(yb(i)-theta(i))/2 ;
   end ;
 /*                    now generate new phi       */
   f = b1new / rangam(seed,a1new) ;
 /*                                               */
 /* new parameters of dist of gamma given others  */
   a2new = a2 + p/2 ;
   b2new = b2 ;
   do i = 1 to p ;
   b2new = b2new + (theta(i)-b)*(theta(i)-b)/2 ;
   end ;
 /*                    now generate new gamma     */
   g = b2new / rangam(seed,a2new) ;
 /*                                               */
 /* new parameters of dist of mu given others     */
   phi0new = 1/( 1/phi0 + p/g ) ;
   beta0new = ( beta0/phi0 + sum( of theta{*} )/g ) * phi0new ;
 /*                    now generate new mu        */
   b = beta0new + sqrt(phi0new)*rannor(seed) ;
 /*                                               */
 /* loop for generating new thetas                */
   do i = 1 to p ;
 /*   first the parameters of dist of theta       */
   gamis = 1/( nn(i)/f + 1/g ) ;
   this = ( nn(i)*yb(i)/f + b/g ) * gamis ;
 /*                    now generate new theta(i)  */
   theta(i) = this + sqrt(gamis)*rannor(seed) ;
   end ;
 /*                                               */
   output ;
   end ;             * this end iteration loop  ;
   cards ;
   6.8  10.3  5.9  6.6  8.5
   2    1     7    3    2
 ;
run ;
  proc print data=a (obs=5) ;
   var f g b t1-t5 ;
   title 'chex131 -- variance components example' ;
   title2 'check the first few observations' ;
run ;
  proc corr data=a noprob cov ;
   var t1 t2 t3 t4 t5 f g b ;
   title2 '' ;
run ;
  proc univariate data=a plot ;
   var f g b ;
run ;
  proc gplot data=a ;
    plot f*g ;
    plot f*b ;
    plot g*b ;
run ;