Example Logs, STAT 705, Fall 2015
===================================

         Part IV, December 2015
         ======================     
                                            12/2/2015

> shoes = data.frame(boy=rep(1:10,2), trt=rep(c("A","B"), 
       rep(10,2)), side=c("L","L","R","L","R","L","L","L","R","L",
       "R","R","L","R","L","R","R","R","L","R"), 
       wear=c(13.2,8.2,10.9,14.3,10.7,6.6,9.5,10.8,8.8,
       13.3,14.0,8.8,11.2,14.2,11.8,6.4,9.8,11.3,9.3,13.6))

> shoes
        boy trt side wear	  ### data from Box, Hunter & Hunter 1978    
      1   1   A    L 13.2    ### discussed on pp.116-118 and p.138
      2   2   A    L  8.2    ### of Venables and Ripley book.
      3   3   A    R 10.9
      4   4   A    L 14.3
      5   5   A    R 10.7
      6   6   A    L  6.6
      7   7   A    L  9.5
      8   8   A    L 10.8
      9   9   A    R  8.8
     10  10   A    L 13.3
     11   1   B    R 14.0
     12   2   B    R  8.8
     13   3   B    L 11.2
     14   4   B    R 14.2
     15   5   B    L 11.8
     16   6   B    R  6.4
     17   7   B    R  9.8
     18   8   B    R 11.3
     19   9   B    L  9.3
     20  10   B    R 13.6

### Permutational illustration for t-test based on centered 
### A-group values and centered B-group values
### if permutations of 1:20 are chosen and restricted so that 
### the R associated with each L receives an opposite symbol.
##(Also enforce that the number of L's assigned A is 2:8.)

> shoes = cbind.data.frame(shoes, wear.ctr = shoes$wear-
      rep(c(mean(shoes$wear[1:10]), mean(shoes$wear[11:20])),
      c(10,10)))

  tmpshoes = shoes[order(shoes$side),]
  inds = 10+match(tmpshoes$boy[1:10],tmpshoes$boy[11:20])
  tmpshoes[11:20,] = tmpshoes[inds,]

> tmpshoes
   boy trt side wear wear.ctr
1    1   A    L 13.2     2.57
2    2   A    L  8.2    -2.43
4    4   A    L 14.3     3.67
6    6   A    L  6.6    -4.03
7    7   A    L  9.5    -1.13
8    8   A    L 10.8     0.17
10  10   A    L 13.3     2.67
13   3   B    L 11.2     0.16
15   5   B    L 11.8     0.76
19   9   B    L  9.3    -1.74
3    1   B    R 14.0     2.96
5    2   B    R  8.8    -2.24
9    4   B    R 14.2     3.16
11   6   B    R  6.4    -4.64
12   7   B    R  9.8    -1.24
14   8   B    R 11.3     0.26
16  10   B    R 13.6     2.56
17   3   A    R 10.9     0.27
18   5   A    R 10.7     0.07
20   9   A    R  8.8    -1.83

tmpshoes[11:20,] = tmpshoes[10+ match(tmpshoes$boy[11:20],      
             tmpshoes$boy[1:10]), ]

> newt1 =  numeric(1000)
  nLvec = rbinom(2000,10,0.5)
  nLvec = nLvec[nLvec %in% 2:8][1:1000]
  for(i in 1:1000) {
      indsA = union(sample(1:10, nLvec[i]),
                 sample(11:20,10-nLvec[i]))
      indsB=setdiff(1:20,indsA)
      newt1[i] = t.test(tmpshoes$wear.ctr[indsA],                   
         tmpshoes$wear.ctr[indsB], var=T)$stat   }

> hist(newt1, nclass=30, prob=T, xlab="t-stat's", main=
   paste("Histogram of t-statistics for Shoe-Data,",
   "\n Permutational based on Residuals"))
  lines(seq(-4,4,length=120), dt(seq(-4,4,length=120),18))

#### Interestingly different !!! Note that the number of
#### allowed configurations is now > sum(choose(10,2:8)) 
### 1002:    VERY limited !!!

## THE POINT OF THIS LAST PICTURE IS THAT THE REMAINING "BUMPS"
## IN THE HISTOGRAM ARE NOW NOT SO MUCH DUE TO SAMPLING AS TO THE 
## GRANULARITY IN THE EXACT PERMUTATIONAL DISTRIBUTION !!

===============================================================
GIBBS SAMPLER INTRODUCTION				12/7/2015
--------------------------

Simplest example:  SIMULATING BIVARIATE NORMAL

Here consider X,Y standard normal variables with correlation rho. Conditionals are 
    f_{Y|X}(y|x) = dnorm(y,rho*x,sqrt(1-rho^2))      and

    f_{X|Y}(x|y) = dnorm(x,rho*y,sqrt(1-rho^2))

Generate RVs X,Y alternately:

   BnormGibbs = function(N, rho) {
       XYmat = array(rnorm(2*N), c(N,2))
       ralt = sqrt(1-rho^2)
       XYmat[1,2] = ralt*XYmat[1,2]+rho*XYmat[1,1]
       for(i in 2:N) {
            XYmat[i,1] = ralt*XYmat[i,1]+ rho*XYmat[i-1,2]            
            XYmat[i,2] = ralt*XYmat[i,2]+ rho*XYmat[i,1] }
       XYmat }

> tmp1 = BnormGibbs(1e4,0.6)
> sqrt(0.09+0.49+2*0.6*0.21)  ## = 0.912
  hist(tmp1[1:300,]%*% c(0.3,0.7), nclass=16, prob=T, ylim=
    c(0,0.6))  curve(dnorm(x,0,0.912),-3,3, add=T, col="red")
         ### agreement not wonderful

> hist(tmp1[301:1e4,]%*% c(0.3,0.7), nclass=50, prob=T, ylim=c(0,0.5))  
  curve(dnorm(x,0,0.912),-4,4, add=T, col="red")
        ### agreement virtually perfect

Second example (from Robert & Casella 1999 book):

      BESAG's AUTOEXPONENTIAL

  f(y1,y2,y3) prop. to.
           exp(-(y1+y2+y3+th1*y1*y2+th2*y1*y3+th3*y2*y3))

    f(y3|y1,y2) = dexp(y3,1+th2*y1+th3*y2)
    f(y2|y1,y3) = dexp(y2,1+th1*y1+th3*y3)
    f(y1|y2,y3) = dexp(y1,1+th1*y2+th2*y3)

    TriAutExp = function(N,th1,th2,th3) {
        Ymat = array(rexp(N*3), c(N,3))
        Y03 = rexp(1)
        Ymat[1,2] = Ymat[1,2]/(1+th1*Ymat[1,1]+th3*Y03)
        Ymat[1,3] = Ymat[1,3]/(1+th2*Ymat[1,1]+th3*Ymat[1,2])
        for(i in 2:N) {
            Ymat[i,1] = Ymat[i,1]/(1+th1*Ymat[i-1,2]+th2*Ymat[i-1,3])            
            Ymat[i,2] = Ymat[i,2]/(1+th1*Ymat[i,1]+th3*Ymat[i-1,3])  
            Ymat[i,3] = Ymat[i,3]/(1+th2*Ymat[i,1]+th3*Ymat[i,2]) }         
        Ymat }

> tmp2 = TriAutExp(1e4,.1,.2,.3)  hist(tmp2[1:500,2], nclass=40,
        col="blue", prob=T)
  plot(hist(tmp2[501:1e4,2], nclass=40, prob=T, plot=F),        
        freq=F, add=T, col="red")
  hist(tmp2[1:3000,2], nclass=40, col="blue", prob=T)
  plot(hist(tmp2[3001:1e4,2], nclass=40, prob=T, plot=F),            
        freq=F, add=T, col="red")
    ### In this example, can try numerical integration to find the
    ###    density at histogram midpoints.

fn1 = function(y) {
          ly = length(y)
          out = numeric(ly)
          for(j in 1:ly)
              out[j] = integrate(function(x,.y)
                   exp(-(x+.y+0.1*x*.y))/(1+0.2*x+0.3*.y),
                   0,Inf,.y=y[j])$val
          out }

> C = 1/integrate(fn1,0,20)$val> C[1] 1.530916

> tmp0 = hist(tmp2[2001:1e4,2], prob=T, xlim=c(0,5), nclass=50, 
           col="cyan")
  curve(C*fn1(x),0,5,add=T)

> qqplot(tmp1[1:1200,]%*% c(0.3,0.7),
         tmp1[1201:1e4,]%*% c(0.3,0.7))
  qqplot(tmp2[1:1000,2],tmp2[1001:1e4,2])

##>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>        12/9/2015                                             

### THIRD EXAMPLE: THE ONE FROM CLASS WEIGHTED-REGRESSION EXAMPLE

Recall
 X[i] ~ Unif[0,1], Y[i] = 3 - 7*x[i] + 
               10*(0.1+abs(x[i]-0.5))^1.5*rnorm(1)

Try prior density   b1,b2 = indep. normal  N(0,10^2)                         
                    tau = 1/sig^2 ~ Gamma(1,0.1)

Then our problem is how to simulate from the posterior 
proportional to

      tau^(n/2) * exp(-(b1^2+b2^2)/2) * 
          exp(-(0.1 + 0.5*sum(w[i]*(y[i]-b1-b2*x[i])^2)*tau)                  

where w[i] = 1/(0.1+abs(x[i]-0.5))^3,  noting that for 
fixed x,w,y data-vectors (which we condition on in taking 
posterior density)

f_{tau|b1,b2} ~ Gamma(1+n/2, 0.1 + 0.5*sum(w[i]*(y[i]-b1-b2*x[i])^2))    
f_{b1,b2|tau} ~ bivariate Normal ((mu1,mu2), Sigma)

with sigma_{12} = rho*sig1*sig2, and 
     tau*sum(w[i])        =  1/(sig1^2*(1-rho^2))'
     tau*sum(w[i]*x[i]^2) =  1/(sig2^2*(1-rho^2))
     tau*sum(w[i]*x[i])   =  -rho/((1-rho^2)*sig1*sig2)
     tau*sum(w[i]*y[i])   =  mu1/(sig1^2*(1-rho^2))
     tau*sum(w[i]*y[i]*x[i]) = mu2/(sig2^2*(1-rho^2)) 

### These equations lead to the MVnorm parameter equations

mu1 =sum(w[i]*y[i])/sum(w[i])
mu2 =sum(w[i]*y[i]*x[i])/sum(w[i])    both of these free of tau
rho = -sum(w[i]*x[i])/sqrt(sum(w[i])*sum(w[i]*x[i]^2))
                                     also free of tau
sig1^2 = { tau*(1-rho^2)*sum(w[i]) }^(-1)
sig2^2 = { tau*(1-rho^2)*sum(w[i]*x[i]^2) }^(-1)

### Simplest expression for (b1,b2) dist is :
## TSum = c(sum(w), sum(w*y), sum(w*y*x), sum(w*x^2), sum(w*x), 
##       sum(w*y^2), length(w))BVNpar = c(Tsum[2]/Tsum[1], 
##       Tsum[3]/Tsum[1], -Tsum[5]/sqrt(Tsum[1]*Tsum[4]),
##       1/(Tsum[1]*(1-BVNpar[3]^2)), 1/(Tsum[4]*(1-BVNpar[3]^2))) 
### these are respectively the expressions for mu1, mu2, rho, 
###      tau*sig1^2, and tau*sig2^2
### Then code Gibbs Sampler alternating between tau and (b1,b2) 
##       as follows:

btauPost = function(N, Tsum, bhat) {
          require(MASS)
          btauArr = array(c(rep(0,N), rnorm(2*N)), c(N,3))
      ### the i'th row will be the i'th round of Gibbs-Sampler
      ### generated  (b1,b2,tau) values
          rho = -Tsum[5]/sqrt(Tsum[1]*Tsum[4])
          muv = Tsum[2:3]/Tsum[1]
          Sig0 = array(c(1/Tsum[1], rep(-Tsum[5]/(Tsum[1]*Tsum[4]),
                2), 1/Tsum[4])/(1-rho^2),c(2,2))
          aux = eigen(Sig0)
          RtSig0 = aux$vec %*% diag(sqrt(aux$val)) %*% t(aux$vec)
      ### symmetric sq rt of Sig0
      ### Initialize using MLE and Infmat
          xmn = Tsum[5]/Tsum[1]
          sxx = Tsum[4]-Tsum[1]*xmn^2
          b2ML = (Tsum[3] - xmn*Tsum[2])/sxx
          b1ML = Tsum[2]/Tsum[1] - b2ML*xmn
          mse = (Tsum[6]-2*b1ML*Tsum[2]-2*b2ML*Tsum[3]+b1ML^2*
              Tsum[1]+2*b1ML*b2ML*Tsum[5]+b2ML^2*Tsum[4])/Tsum[7]             
          Infmat = array(c(Tsum[1],Tsum[5],0,Tsum[5],Tsum[4],
              0,0,0,0.5*Tsum[7]/mse),c(3,3))/mse
          Vmat = solve(Infmat)
          btauArr[1,1:2] = mvrnorm(1,c(b1ML,b2ML),Vmat[1:2,1:2])     
          btauArr[1,3] = 1/(mse*exp(rnorm(1)*sqrt(Vmat[3,3])))   
          for(i in 2:N) {
             btauArr[i,1:2] = c(b1ML,b2ML) + (RtSig0 %*% 
                 btauArr[i,1:2])/sqrt(btauArr[i-1,3])
             btauArr[i,3] = rgamma(1,1+Tsum[7]/2,0.1+0.5*            
                (Tsum[6]-2*btauArr[i,1]*Tsum[2]-2*btauArr[i,2]*
                 Tsum[3]+btauArr[i,1]^2*Tsum[1]+2*btauArr[i,1]*
                 btauArr[i,2]*Tsum[5]+ btauArr[i,2]^2*Tsum[4])) }
          btauArr }

  TSum0 = c(sum(Dset$wt[1:100]), sum((Dset$wt*
              Dset$y)[1:100]), sum((Dset$wt*Dset$y*Dset$x)[1:100]), 
              sum((Dset$wt*Dset$x^2)[1:100]), 
              sum((Dset$wt*Dset$x)[1:100]),
              sum((Dset$wt*Dset$y^2)[1:100]), 100)

  tmp = btauPost(1e4, TSum0, tmp1$coef)  ### seems OK

  unix.time({ tmp = btauPost(1e5, TSum0, tmp1$coef) })
   user  system elapsed
    1.98    0.00    1.98 

### Now we examine the three histograms

  hist(tmp[,1], nclass=50, xlab="b1 parameter", prob=T,
    main="Gibbs-Sampler Posterior Density, b1", ylim=c(0,1.5))
    abline(v=tmp1$coef[1], col="red")
  curve(dnorm(x,tmp1$coef[1], sqrt(solve(Infmat)[1,1])),       
      add=T, col="blue") 

### We can see clearly that the "asymptotic" form of the 
### posterior is not quite accurate in this setting with n=100, 
### and the same is true for the b2 posterior.
### Pictures are saved as Gibbs_b1.pdf, Gibbs_b2.pdf

  hist(tmp[,2], nclass=50, xlab="b2 parameter", prob=T,
    main="Gibbs-Sampler Posterior Density, b2", ylim=c(0,0.7))
  abline(v=tmp1$coef[2], col="red")
  curve(dnorm(x,tmp1$coef[2],sqrt(solve(Infmat)[2,2])),       
        add=T, col="blue") 

### Finally, the histogram and "theoretical" density for sig^2

>  hist(1/tmp[,3], nclass=200, xlab="sig^2 parameter", xlim=
       c(45,155), ylim=c(0,.04), prob=T, main=
       "Gibbs-Sampler Posterior Density, sig^2")
   abline(v=summary(tmp1)$sig^2*98/100, col="red")
   curve(dnorm(x,summary(tmp1)$sig^2*98/100,
         sqrt(solve(Infmat)[3,3])), add=T, col="blue")