#  Gibbs sampling estimation of variances
#  Same model

# set up arrays to accumulate sample variances
CGV = c(0)   # for contemporary groups
AV = CGV     # for additive genetic
RV = CGV     # for residual
VRcg = CGV   # for ratio of residual to contemporary group variances
VRadd = CGV  # for ratio of residual to additive variances
VRh2 = CGV   # for heritability
#  Alterations made to MME include Gibbs Sampling
#  Just 100 samples will be computed here

#  Begin Samples
resid=700   # starting value of residual variance
iter = 0
while(iter < 10000){
# New Contemporary group solutions (Random Factor in Model)
#  Must Add variance ratio, solutions should sum to zero
iter = iter + 1
rhs = y - sage[agefac] - (bhat * cov) - sanm[anfac]
rhs
scg = tapply(rhs,cgfac,sum)
dd = as.numeric(tapply(rhs,cgfac,length))
dd = dd + ratcg
dd
n = length(scg)  # number of contemporary groups
v = rnorm(n,0,1)*sqrt(resid/dd)
scg = (scg/dd) + v
scg

#  Estimate CG variance, accumulate in a vector
S = as.numeric(t(scg) %*% scg)
xd = rchisq(1,n)
cgvar = S/xd
CGV = c(CGV, cgvar)

# New Age Factor solutions (Fixed factor in Model)
rhs = y - scg[cgfac] - (bhat*cov) - sanm[anfac]
sage = tapply(rhs,agefac,mean)
dd = as.numeric(tapply(rhs,agefac,length))
n = length(sage)
v = rnorm(n,0,1)*sqrt(resid/dd)
sage = sage + v

# New marker regression coefficient
rhs = y - scg[cgfac] - sage[agefac] - sanm[anfac]
bhat = xx %*% covt %*% as.matrix(rhs)
v = rnorm(1)*sqrt(resid*xx)
bhat = bhat + v

# New animal genetic solutions (random factor in model
#  Must account for relationships among animals
rhs = y - scg[cgfac] - sage[agefac] - (bhat*cov)
rhs
arhs = rep(0,nanim+1)
k2=b2$aaid
arhs[k2]=rhs

ij = kord[1]
kprev = aaa[ij]
for(j in 1:nped) {     # Use coded pedigree list
   i = kord[j]
   ka = aaa[i]
   ks = sss[i]
   kd = ddd[i]
  if(ka != kprev){
   d = bii[kprev]*ratan
   x = zz[kprev] + ratan*adiag[kprev]
   anew = arhs[kprev]/x
   v = rnorm(1)*sqrt(resid/x)
   sanm[kprev] = anew + v
   }
   kprev = ka
   if(cods[i]==0){
         d = bii[ka]*ratan
         arhs[ka]=arhs[ka] + 0.5*d*(sanm[ks]+sanm[kd])
        } else    # It is important not to put 'else' by itself
        { d = bii[ks]*ratan
         arhs[ka] = arhs[ka] + d*(0.5*sanm[ks]-0.25*sanm[kd]) }
 }
#  Take care of last animal
   d = bii[kprev]*ratan
   x = zz[kprev] + ratan*adiag[kprev]
   anew = arhs[kprev]/x
   v = rnorm(1)*sqrt(resid/x)
   sanm[kprev] = anew + v
sanm

S = 0
for( i in 1:nanim){
    MS = sanm[i] - 0.5*(sanm[ssid[i]] + sanm[ddid[i]])
    S = S + MS*MS*bii[i]
   }
xd = rchisq(1,nanim)
addvar = S/xd
AV = c(AV, addvar)

#  Estimate residual variance
rhs = y - scg[cgfac] - sage[agefac] - (bhat*cov) - sanm[anfac]
S = as.numeric(t(rhs) %*% rhs)
n = length(rhs)
xd = rchisq(1,n)
resid = S/xd
RV = c(RV, resid) 
ratcg = resid/cgvar
if(ratcg > 25){ratcg = pricg}
ratan = resid/addvar
if(ratan > 25){ratan = prian}
 
h2 = addvar/(addvar + resid + cgvar)
VRcg = c(VRcg, ratcg)
VRadd = c(VRadd, addvar)
VRh2 = c(VRh2, h2)   }