R Log for Solution of HW 13, Stat 705
=====================================

## Simulate 10 x 40 array of epsilon ~ N(0,1)
## 40-vec U ~ N(0,2), 10-vec alpha ~ Unif(8,12)

eps = array(rnorm(400), c(10,40))
Uvec = sqrt(2)*rnorm(40)
alpha = runif(10,8,12)
Xarr = outer(alpha,Uvec,"+") + eps
ahat = c(Xarr %*% rep(1,40))/40
SSW = sum((-ahat + Xarr)^2)
SSBC = 10*39*var(c(t(Xarr) %*% rep(.1,10)))
sgsqhat = (SSW-SSBC)/(40*9)
sgsquhat = (SSBC/40-sgsqhat)/10
> round(rbind(tru=alpha,est=ahat),3)
      [,1]   [,2]   [,3]  [,4]  [,5]  [,6]   [,7]   [,8]  [,9]  [,10]
tru 11.635 10.786 10.982 9.293 8.768 8.228 11.437 11.498 8.777 11.279
est 11.642 10.917 10.875 9.303 8.997 8.212 11.180 11.618 9.085 11.092
> c(sgsqhat, sgsquhat)
[1] 0.8608038 2.6643837    ### compare with 1, 2

EMitrN = function(sgsq,sgsqu,B,M) {
      gam = sgsqu/(sgsqu+sgsq/B)
      sgsqun = (gam/B)*(sgsq+SSBC*gam/M)
      sgsqn = sgsqun + (SSW-2*gam*SSBC)/(M*B)
      c(sgsqn, sgsqun) }

> EMitrN(.1,.1,10,40)
[1] 0.8065455 2.2822017

## Start with initial value th = c(ahat,1e-2,1e-2)

parm = rep(1e-2,2)
parm=EMitrN(parm[1],parm[2],10,40)  ## =  0.7983636 2.2740199
parm=EMitrN(parm[1],parm[2],10,40)  ## =  0.855016 2.644180

for (i in 1:10) {
   parm=EMitrN(parm[1],parm[2],10,40)  
   cat(round(parm,6),"\n") }
0.860245 2.663681 
0.860749 2.664393 
0.860798 2.664389 
0.860803 2.664385 
0.860804 2.664384 
0.860804 2.664384 
0.860804 2.664384 
0.860804 2.664384 
0.860804 2.664384 
0.860804 2.664384   ### Converged very quickly !!


### Now do it all 100 times, both with nlm and EM:

Minusllk = function(parm, ssw, ssbc, B=10, M=40) {
        varsum = B*parm[2]+parm[1]
        gam = B*parm[2]/varsum
        (M*(B-1)/2)*log(parm[1])+(M/2)*log(varsum)+
                (ssw-gam*ssbc)/(2*parm[1])   
     }

parm = array(0, c(100,2))
llkseq  = numeric(10)
for(i in 1:100) {
        Uvec = sqrt(2)*rnorm(40)
        eps = array(rnorm(400), c(10,40))
        alpha = runif(10,8,12)
        Xarr = outer(alpha,Uvec,"+") + eps
        ahat = c(Xarr %*% rep(1,40))/40
        SSW = sum((-ahat + Xarr)^2)
        SSBC = 10*39*var(c(t(Xarr) %*% rep(.1,10)))
        sgsqhat = (SSW-SSBC)/(40*9)
        sgsquhat = (SSBC/40-sgsqhat)/10
        parm[i,]=rep(1e-2,2)
        llkseq[1] = -Minusllk(parm[i,],SSW,SSBC)
        for(j in 2:10) {  parm[i,] = 
             EMitrN(parm[i,1],parm[i,2],10,40)
             llkseq[j] = -Minusllk(parm[i,],SSW,SSBC) }
        if(min(diff(llkseq)) < -1e-12) cat("!Decr logLik!! \n")
        parm0 = nlm(Minusllk, parm[i,], ssw=SSW,ssbc=SSBC)$est
        if (sum(abs(c(sgsqhat,sgsquhat)-parm[i,])) > 1e-4 |
             sum(abs(c(sgsqhat,sgsquhat)-parm0)) > 1e-4 ) 
               cat("!Unequal Est's!! \n")
        }

### No printing occurred.

> summary(parm)
       V1               V2
 Min.   :0.8283   Min.   :0.8754
 1st Qu.:0.9207   1st Qu.:1.5376
 Median :0.9650   Median :1.8957
 Mean   :0.9707   Mean   :1.8943
 3rd Qu.:1.0083   3rd Qu.:2.1838
 Max.   :1.1880   Max.   :3.5297
### Occasional poor estimates are OK !!
###   Note also that the MLE's of variances are 
###     somewhat biased below in this example !!