rmvnorm<-function(n,mu,Sigma,chol.Sigma=chol(Sigma)) {
  E<-matrix(rnorm(n*length(mu)),n,length(mu))
  t(  t(E%*%chol.Sigma) +c(mu))
                              }


rwish<-function(S0,nu){       # sample from a Wishart distribution
  sS0<-chol(S0)
  Z<-matrix(rnorm(nu*dim(S0)[1]),nu,dim(S0)[1])%*%sS0
  t(Z)%*%Z
                       }




rF.fc<-function(Z,G,E,C,s2)
{ 
  R<-dim(G)[2] ; m<-dim(Z)[1] 
  v.F<- solve( solve(C) +  t(G)%*%G/s2 )
  m.F<- v.F%*%( c(solve(C)%*%c(E)) + t(G)%*%t(Z)/s2 )
  X<-matrix(rnorm(m*R),m,R)
  t( t(X%*%chol(v.F))  + m.F ) 
}



rC.fc<-function(U,mu)
{
  R<-dim(U)[2] ; m<-dim(U)[1]
  E<-t( t(U)-c(mu))
  solve(rwish(solve( diag(nrow=R)*sqrt(R) +t(E)%*%E ),m+R+2))
}


XB<-function(X, B) 
{
    tmp <- matrix(0, nrow = dim(X)[1], ncol = dim(X)[2])
    for (k in seq(1, length(B), length = length(B))) {
        tmp <- tmp + B[k] * X[, , k]
    }
    tmp
}


rB.fc<-function(Z,xt,ixtx,s2)
{
  vr<-ixtx/( 1/prod(dim(Z))+1/s2) 
  mn<-vr%*%( xt%*%c(Z) )
  rmvnorm(1,mn,vr)
}


rE.fc<-function(U,C)
{
  R<-dim(U)[2] ;  m<-dim(U)[1]
  v.mu<-solve( diag(1/sqrt(m),nrow=R) + m*solve(C) )
  m.mu<-v.mu%*%( solve(C)%*%apply(U,2,mean) )
  rmvnorm(1, m.mu,v.mu)
}

rZ.fc<-function(EZ,Y)
{
  for(y in sort(unique(c(Y))) )
  {
    ir<-(Y==y & !is.na(Y))
    lb<-suppressWarnings( max(Z[Y<y],na.rm=TRUE) )
    ub<-suppressWarnings( min(Z[Y>y],na.rm=TRUE) )
    z<-qnorm(runif(sum(ir),pnorm(lb, EZ[ir]), pnorm(ub,EZ[ir],1)),EZ[ir])
    z[z == Inf] <- lb
    z[z == -Inf] <- ub
    Z[ir]<- z
  } 
  ir<-is.na(Y)
  Z[ir] <- rnorm(sum(ir), EZ[ir], 1)

  Z
}