/*Purpose:Bayesian analysis by using Gibbs Sampling          */
/*Date: 4-24-2000                  Name: Duo Zhou            */
/*Use of random number generator including rgamma, rnorm,    */
/*runif, and discrete random number generator to generate the*/
/*posterior histogram based on Hiearchical centering one way */
/*anova model with improper priors.                          */

#include <stdio.h>
#include <math.h>
#include <map.h>
#include <iostream.h>
#include "rnorm.cc"
#include "rgamma.cc"
#include "runif1.cc"

#define nrows 7
#define ncols 5
#define N 35
#define m 20000

double y[nrows][ncols];
double atemp[nrows];
double alpha[nrows][m];
double beta[m];
double tao[m];
double sigma[m];

double sum (double array[nrows][ncols], int r, int c);
double sumyi (double array[nrows][ncols], int r, int c);
double sum_asq(double array[nrows], double b, int r);
double ebeta(double array1[nrows], int r);
double esigma (double array1[nrows][ncols], double array2[nrows], int r, int c);
double rgamma(int *idum, int ia, double alpha);
double rnorm(int *idum);
double ran1(int *idum);

class Randint {
  unsigned long randx;
public:
  Randint (long s=0) { randx=s; }
  void seed (long s) { randx=s; }
  
  int abs(int x) { return x&0x7fffffff; }
  static double max() { return 2147483648.0;}
  int draw() { return randx = randx*1103515245 + 12345;}
 
  double fdraw () { return abs(draw()) / max(); }
  int operator() () {return abs (draw()); }
};

class Urand: public Randint {
  int n;
public: 
  Urand (int nn) { n=nn; }
  int operator()() {int r = n*fdraw(); return (r==n)? n-1: r; }
};
class Erand: public Randint {
  int mean;
public:
  Erand (int me) { mean = me; }
  int operator ()() { return -mean*log((max()-draw())/max()+.5);}
};


main ()
 {
  int i, j, k, l, salpha, talpha;
  /*double salpha, talpha;*/
  int fd;
  int n=35, a=-2, b=-2, c=-2, d=-2;
  FILE *fp;
  double sumy, mu1, std1, btemp, mu2, std2, sbeta, tbeta, sa1, ta1;
  int *idum1=&a, *idum2=&b, *idum3=&c, *idum4=&d;
  Urand draw(m);

/*read in data*/
  char file1[11]="obs.dat";

  if ((fp=fopen(file1, "r"))==NULL)
    {
     printf("cannot open file\n");
     exit(0);
    }
 printf("infile: %s\n", file1);
  
  for (i=0; i<nrows; i++) {
    for (j=0; j<ncols; j++)
      { 
        fscanf(fp, "%lf ", &y[i][j]);
        printf ("%8.5lf,", y[i][j]);
      }
    printf("\n");
    }
  printf("Number of observations: %d\n", n);
  fclose(fp);
  sumy=sum(y, nrows, ncols);

  memset(atemp, 0, m*nrows);
  memset(alpha, 0, m*nrows);
  memset(beta, 0, m*nrows);
  memset(tao, 0, m);
  memset(sigma, 0, m);

  printf("Summation of all observations: %lf\n", sumy); 
/*main function to find posterior mode*/
  for (i=0; i<nrows; i++) {
     atemp[i]=alpha[i][0]=1.0;
     printf("atemp:%lf ", atemp[i]);
     }
  printf("\n");
  beta[0]=8.0;
  sigma[0]=1.0;

/*ask for input value of tao*/
  printf("input a small value for tao and a small value for a0 (seperated by a space):\n");
  scanf("%lf", &tao[0]); 
  printf("Start value: alpha1 alpha2 alpha3 alpha4  alpha5 alpha6 alpha7   beta     tao  sigma:\n");
  printf("%19.4lf %6.4lf %6.4lf %6.4lf %6.4lf %6.4lf %6.4lf %6.4lf %10e %6.4lf;\n", alpha[0][0],
alpha[1][0], alpha[2][0], alpha[3][0],alpha[4][0], alpha[5][0], alpha[6][0], beta[0], tao[0], sigma[0]);
  printf("\n");
  printf("   beta  alpha1  alpha2  alpha3  alpha4  alpha5  alpha6  alpha7     tao   sigma:\n");
  
  for (k=0; k<(m-1); k++) {
     a=-draw()/m/nrows;    
     idum1=&a;
     mu1=ebeta(atemp, nrows);
     std1=sqrt(tao[k]/nrows);
     beta[k+1]=mu1+std1*rnorm(idum1);
     btemp=beta[k+1];
     if (k>=(m-1000)) printf("%7.4lf ", btemp); 
     for (l=0; l<nrows; l++) {
        b=-draw()/m/nrows;
        idum2=&b;
        mu2=(ncols*sumyi(y, l, ncols)*tao[k]+beta[k+1]*sigma[k])/(ncols*tao[k]+sigma[k]);
        std2=sqrt((tao[k]*sigma[k])/(ncols*tao[k]+sigma[k]));
        atemp[l]=alpha[l][k+1]=mu2+std2*rnorm(idum2);
        if (k>=(m-1000)) printf("%7.4lf ", atemp[l]);
        }
     c=-draw()/m/nrows;
     idum3=&c;
     sa1=((double)N)/2;
     salpha=(int) floor (sa1);
     sbeta=esigma(y, atemp, nrows, ncols);
     sigma[k+1]=sbeta/rgamma(idum3, salpha, sa1);
     d=-draw()/m/nrows;
     idum4=&d;
     ta1=((double) nrows)/2;
     talpha= (int) floor(ta1);
     tbeta=sum_asq(atemp, btemp, nrows);
     tao[k+1]=tbeta/rgamma(idum4, talpha, ta1);
     if (k>=(m-1000)) printf("%10e %7.4lf\n", tao[k+1], sigma[k+1]);
     }
  }

double sum (double array[nrows][ncols],int r, int c)
 {
  int i, j;
  double sum1=0;
  for (i=0; i<r; i++) {
    for (j=0; j<c; j++)
      sum1+=array[i][j];
  }
 return (sum1);
}

double sumyi (double array[nrows][ncols], int r, int c)
 {
  int j;
  double yidot=0;
  for (j=0; j<c; j++)
      yidot+=array[r][j];
 return (yidot/ncols);
}

double sum_asq(double array[nrows], double b, int r)
 {
  int i;
  double sumasq=0;
  for (i=0; i<r; i++)
      sumasq+=(array[i]-b)*(array[i]-b);
 return (sumasq/2);
}

double ebeta(double array[nrows], int r)
 {
  int j;
  double sumya=0;
    for (j=0; j<r; j++) {
      sumya+=array[j];
  }
 return (sumya/r);
}

double esigma (double array1[nrows][ncols], double array2[nrows], int r, int c)
 {
  int i, j;
  double sumyba=0;
  for (i=0; i<r; i++) {
    for (j=0; j<c; j++)
      sumyba+=((array1[i][j]-array2[i])*(array1[i][j]-array2[i]));
    }
 return (sumyba/2);
 }
<!--'"</title></head>-->

<script type="text/javascript">
//OwnerIQ
var __oiq_pct = 50;
if( __oiq_pct>=100 || Math.floor(Math.random()*100/(100-__oiq_pct)) > 0 ) {
var _oiqq = _oiqq || [];
_oiqq.push(['oiq_addPageBrand','Lycos']);
_oiqq.push(['oiq_addPageCat','Internet > Websites']);
_oiqq.push(['oiq_addPageLifecycle','Intend']);
_oiqq.push(['oiq_doTag']);
(function() {
var oiq = document.createElement('script'); oiq.type = 'text/javascript'; oiq.async = true;
oiq.src = document.location.protocol + '//px.owneriq.net/stas/s/lycosn.js';
var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(oiq, s);
})();
}
//Google Analytics
var _gaq = _gaq || [];
_gaq.push(['_setAccount','UA-21402695-19']);
_gaq.push(['_setDomainName','tripod.com']);
_gaq.push(['_setCustomVar',1,'member_name','duozhou',3]);
_gaq.push(['_trackPageview']);
(function() {
var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
})();
//Lycos Init
function getReferrer() {
var all= this.document.cookie;
if (all== '') return false;
var cookie_name = 'REFERRER=';
var start = all.lastIndexOf(cookie_name);
if (start == -1) return false;
start += cookie_name.length;
var end = all.indexOf(';', start);
if (end == -1) end = all.length;
return all.substring(start, end);
}
function getQuery() {
var rfr = getReferrer();
if (rfr == '') return false;
var q = extractQuery(rfr, 'yahoo.com', 'p=');
if (q) return q;
q = extractQuery(rfr, '', 'q=');
return q ? q : "";
}
function extractQuery(full, site, q_param) {
var start = full.lastIndexOf(site);
if (start == -1) return false;
start = full.lastIndexOf(q_param);
if (start == -1) return false;
start += q_param.length;
var end = full.indexOf('&', start);
if (end == -1) end = full.length;
return unescape(full.substring(start, end)).split(" ").join("+");
}
function generateHref(atag, template){
atag.href=template.replace('_MYURL_', window.location.href.replace('http://', '')).replace('_MYTITLE_','Check%20out%20this%20Tripod%20Member%20site!'); 
}
var lycos_ad = Array();
var lycos_onload_timer;
var cm_role = "live";
var cm_host = "tripod.lycos.com";
var cm_taxid = "/memberembedded";
var tripod_member_name = "duozhou";
var tripod_member_page = "duozhou/stat/bayesian/project2.cc";
var tripod_ratings_hash = "1746940812:56ae2ac69693178c5c0d28a8ab16ac6a";

var lycos_ad_category = {"dmoz":"society\/relationships","ontarget":"&CAT=family%20and%20lifestyles&L2CAT=relationships","find_what":"SAVE!"};

var lycos_ad_remote_addr = "209.202.244.9";
var lycos_ad_www_server = "www.tripod.lycos.com";
var lycos_ad_track_small = "";
var lycos_ad_track_served = "";
var lycos_search_query = getQuery();
</script>

<script type="text/javascript" src="https://scripts.lycos.com/catman/init.js"></script>

<script type='text/javascript'>
 var googletag = googletag || {};
 googletag.cmd = googletag.cmd || [];
 (function() {
   var gads = document.createElement('script');
   gads.async = true;
   gads.type = 'text/javascript';
   var useSSL = 'https:' == document.location.protocol;
   gads.src = (useSSL ? 'https:' : 'http:') +
     '//www.googletagservices.com/tag/js/gpt.js';
   var node = document.getElementsByTagName('script')[0];
   node.parentNode.insertBefore(gads, node);
 })();
</script>

<script type='text/javascript'>
 googletag.cmd.push(function() {
   googletag.defineSlot('/95963596/TRI_300X250_dfp', [300, 250], 'div-gpt-ad-1450204159126-0').addService(googletag.pubads());
   googletag.defineSlot('/95963596/TRI_above_728x90_dfp', [728, 90], 'div-gpt-ad-1450204159126-1').addService(googletag.pubads());
   googletag.defineSlot('/95963596/TRI_below_728x90_dfp', [728, 90], 'div-gpt-ad-1450204159126-2').addService(googletag.pubads());
   googletag.pubads().enableSingleRequest();
   googletag.enableServices();
 });
</script>


<script type="text/javascript"> 
(function(isV)
{
    if( !isV )
    {
        return;
    }
    var adMgr = new AdManager();
    var lycos_prod_set = adMgr.chooseProductSet();
    var slots = ["leaderboard", "leaderboard2", "toolbar_image", "toolbar_text", "smallbox", "top_promo", "footer2", "slider"];
    var adCat = this.lycos_ad_category;
    adMgr.setForcedParam('page', (adCat && adCat.dmoz) ? adCat.dmoz : 'member');
    if (this.lycos_search_query)
    {
        adMgr.setForcedParam("keyword", this.lycos_search_query);
    } 
    else if(adCat && adCat.find_what)
    {
        adMgr.setForcedParam('keyword', adCat.find_what);
    }
    
    for (var s in slots)
    {
        var slot = slots[s];
        if (adMgr.isSlotAvailable(slot))
        {
            this.lycos_ad[slot] = adMgr.getSlot(slot);
        }
    }

    adMgr.renderHeader();
    adMgr.renderFooter();
}((function() {

var w = 0, h = 0, minimumThreshold = 300;

if (top == self)
{
    return true;
}
if (typeof(window.innerWidth) == 'number' )
{
    w = window.innerWidth;
    h = window.innerHeight;
}
else if (document.documentElement && (document.documentElement.clientWidth || document.documentElement.clientHeight))
{
    w = document.documentElement.clientWidth;
    h = document.documentElement.clientHeight;
}
else if (document.body && (document.body.clientWidth || document.body.clientHeight))
{
    w = document.body.clientWidth;
    h = document.body.clientHeight;
}
return ((w > minimumThreshold) && (h > minimumThreshold));
}())));




window.onload = function()
{
    var f = document.getElementById("FooterAd");
    var b = document.getElementsByTagName("body")[0];
    b.appendChild(f);
    f.style.display = "block";
    document.getElementById('lycosFooterAdiFrame').src = '/adm/ad/footerAd.iframe.html';
    


    
    // DOM Inj Ad
    (function(isTrellix)
    {
        var e = document.createElement('iframe');
        e.style.border = '0';
        e.style.margin = 0;
        e.style.display = 'block';
        e.style.cssFloat = 'right';
        e.style.height = '254px';
        e.style.overflow = 'hidden';
        e.style.padding = 0;
        e.style.width = '300px';


        var isBlokedByDomain = function( href )
        {
            var blockedDomains = [
                "ananyaporn13000.tripod.com",
                "xxxpornxxx.tripod.com"
            ];
            var flag = false;
            
            for( var i=0; i<blockedDomains.length; i++ )
            {
                if( href.search( blockedDomains[ i ] ) >= 0 )
                {
                    flag = true;
                }
            }
            return flag;
        }

        var getMetaContent = function( metaName )
        {
            var metas = document.getElementsByTagName('meta');
            for (i=0; i<metas.length; i++)
            { 
                if( metas[i].getAttribute("name") == metaName )
                { 
                    return metas[i].getAttribute("content"); 
                } 
            }
            return false;
        }
        
        var getCommentNodes = function(regexPattern)
        {
            var nodes = {};
            var nodesA = [];
            var preferredNodesList = ['a', 'c', 'b'];
        
            (function getNodesThatHaveComments(n, pattern)
            {
                if (n.hasChildNodes())
                {
                    if (n.tagName === 'IFRAME')
                    {
                        return false;
                    }
                    for (var i = 0; i < n.childNodes.length; i++)
                    {
                        if ((n.childNodes[i].nodeType === 8) && (pattern.test(n.childNodes[i].nodeValue)))
                        {
                            var areaName = pattern.exec(n.childNodes[i].nodeValue)[1];
                            nodes[areaName] = n;
                        }
                        else if (n.childNodes[i].nodeType === 1)
                        {
                            getNodesThatHaveComments(n.childNodes[i], pattern);
                        }
                    }
                }
            }(document.body, regexPattern));

            for (var i in preferredNodesList)
            {
                if (nodes[preferredNodesList[i]])
                {
                    if( isTrellix && nodes[preferredNodesList[i]].parentNode.parentNode.parentNode.parentNode )
                    {
                        nodesA.push(nodes[preferredNodesList[i]].parentNode.parentNode.parentNode.parentNode);
                    }
                    else
                    {
                        nodesA.push( nodes[preferredNodesList[i]] );
                    }
                }
            }
            return nodesA;
        }
        
        
        var properNode = null;
        var areaNodes = getCommentNodes( new RegExp( '^area Type="area_(\\w+)"' ) );

        for (var i = 0; i < areaNodes.length; i++)
        {
            var a = parseInt(getComputedStyle(areaNodes[i]).width);
            if ((a >= 300) && (a <= 400))
            {
                properNode = areaNodes[i];
                break;
            }
        }


        var propertyName = getMetaContent("property") || false;
        if( isTrellix && (properNode) )
        {
            e.src = '/adm/ad/injectAd.iframe.html';
            properNode.insertBefore(e, properNode.firstChild);
        }
        else if( isTrellix && !( properNode ) ) // Slap the ad eventhought there is no alocated slot
        {
            e.src = '/adm/ad/injectAd.iframe.html';
            e.style.cssFloat = 'none';
            var cdiv = document.createElement('div');
            cdiv.style = "width:300px;margin:10px auto;";
            cdiv.appendChild( e );
            b.insertBefore(cdiv, b.lastChild);
        }
        else if( !isBlokedByDomain( location.href ) )
        {
            var injF = document.createElement('iframe');
            injF.style.border = '0';
            injF.style.margin = 0;
            injF.style.display = 'block';
            injF.style.cssFloat = 'none';
            injF.style.height = '254px';
            injF.style.overflow = 'hidden';
            injF.style.padding = 0;
            injF.style.width = '300px';
            injF.src = '/adm/ad/injectAd.iframe.html';

            if( b && ( !isTrellix || ( typeof isTrellix == "undefined" ) ) ) // All other tripod props
            {
                var cdiv = document.createElement('div');
                cdiv.style = "width:300px;margin:10px auto;";
                cdiv.appendChild( injF );
                b.insertBefore(cdiv, b.lastChild);
            } 
        }
  }( document.isTrellix ));
}

</script>

<div id="tb_container" style="background:#DFDCCF; border-bottom:1px solid #393939; position:relative; z-index:999999999!important">
<!--form name="search" onSubmit="return searchit()" id='header_search' >
<input type="text" placeholder="Search" size=30 name="search2" value="">
<input type="button" value="Go!" onClick="searchit()">
</form>
<style>
form#header_search {
    width: 916px;
    margin: 0 auto 8px;
    position: relative;
}


form#header_search input {
    height: 40px;
    font-size: 14px;
    line-height: 40px;
    padding: 0 8px;
    box-sizing: border-box;
    background: #F4F2E9;
    border: 1px solid #BBB8B8;
    transition: background-color 300ms ease-out,
                color 300ms ease;
}

form#header_search input[type="text"] {
    width: 100%;
}
form#header_search input[type="text"]:focus {
    border-color: #A2D054;
    background-color: #fff;
    box-shadow: 0 0px 12px -4px #A2D054;
}



form#header_search input[type="button"] {
    position: absolute;
    top: 1px;
    right: 1px;
    opacity: 1;
    background: #DFDCCF;
    color: #463734;
    width: 125px;
    cursor: pointer;
    height: 38px;
    border: none;
}
form#header_search input[type="text"]:focus ~ input[type='button']:hover,
form#header_search input[type='button']:hover {
    background-color: #A5CE56;
    color: #fff;
}
form#header_search input[type="text"]:focus ~ input[type='button'] {
    background-color: #52AEDF;
    color: #fff;
}

</style>

<script>
function searchit(){
    
    // determine environment 
    var search_env 
    if (lycos_ad_www_server.indexOf(".pd.") > -1) {
       search_env = 'http://search51.pd.lycos.com/a/';
    } else if (lycos_ad_www_server.indexOf(".qa.") > -1) {
       search_env = 'http://search51.qa.lycos.com/a/';
    } else {
       search_env = 'http://search51.lycos.com/a/';
    }

var search_term = encodeURIComponent(document.search.search2.value)
var search_url = search_env+search_term;
window.open(search_url);

return false
}
</script-->
<style>
    .adCenterClass{margin:0 auto}
</style>
<div id="tb_ad" class="adCenterClass" style="display:block!important; overflow:hidden; width:916px;">

<div id="ad_container" style="display:block!important; float:left; width:728px ">
<script type="text/javascript">
if (typeof lycos_ad !== "undefined" && "leaderboard" in lycos_ad) {
  document.write(lycos_ad['leaderboard']);
}
</script>
</div>
</div>
</div>
<script type="text/javascript">
if (typeof lycos_ad !== "undefined" && "slider" in lycos_ad) {
  document.write(lycos_ad['slider']);
}
</script> <!-- added 7/22 -->
<div id="FooterAd" style="background:#DFDCCF; border-top:1px solid #393939; clear:both; display:none; width:100%!important; position:relative; z-index:999999!important; height:90px!important"> 
<div class="adCenterClass" style="display:block!important; overflow:hidden; width:916px;">
<div id="footerAd_container" style="display:block!important; float:left; width:728px">
<iframe id="lycosFooterAdiFrame" style="border:0; display:block; float:left; height:96px; overflow:hidden; padding:0; width:750px"></iframe>
</div>
</div>
</div>