Jeremy' s IQ Example .  Hierarchical Normal Model

data = {102,  112,  
        96,  109,  98}

Mean[data] // N

s2f = 80
beta = 110
tau2  =  120
p = 5
s2 = 
      Apply[Plus,  (data - Mean[data])^2]
mstar  =  
      data  -  
        s2f/(s2f + s2p) (data - 
              Mean[data]) -   (s2f + s2p)/(p tau2  +  s2f + 
                s2p) (Mean[data] - beta) {1, 1, 1, 1, 1}

con = 
      NIntegrate[
        Sqrt[tau2] Exp[-1/2 ( 
                  s2/(s2f  +  s2p) +  
                    p (Mean[data] - beta)^2/(p tau2  +  s2f  + 
                            s2p))]/((s2f + s2p)^((p - 1)/
                      2) (p tau2  +  s2f  +  s2p)^(1/
                      2)), {s2p, 0, Infinity}]

den  =  
      Sqrt[tau2] Exp[-1/2 ( 
                s2/(s2f  +  s2p) +  
                  p (Mean[data] - beta)^2/(p tau2  +  s2f  + 
                          s2p))]/(0.00227811 (s2f + 
                    s2p)^((p - 1)/2) (p tau2  +  s2f  +  
                    s2p)^(1/2))
Plot[den, {s2p, 0, 1000}]

aaa = mstar den

mup  =  
      Table[NIntegrate[aaa[[i]], {s2p, 0, Infinity}], {i, 5}]

I5  =  
        IdentityMatrix[5];

ones = Table[1, {5}];
J5  =   Transpose[{ones}]  .  {ones}

Vp  =  
        s2f  I5  -  
          s2f^2/(s2f  +  
                  s2p)  I5   +   s2f^2 tau2 /((s2f  +  
                      s2p) (s2f  +  s2p + 
                      p tau2) ) J5  +  
 
          Transpose[{mstar  -  mup}]  .   {mstar  -  
                mup} ;

vvv  =  den Vp;

Vpost  =  
      Table[NIntegrate[vvv[[i, j]], {s2p, 0, Infinity}], {i, 5}, {j, 
          5}]