################################################################### # # R source code providing generalized bilinear mixed effects modeling # for social network data. For more information on the model fitting # procedure, see # "Bilinear Mixed-Effects Models for Dyadic Data" (Hoff 2003) # http://www.stat.washington.edu/www/research/reports/2003/tr430.pdf # # This research was supported by ONR grant N00014-02-1-1011 # # Copyright (C) 2003 Peter D. Hoff # # This program is free software; you can redistribute it and/or # modify it under the terms of the GNU General Public License # as published by the Free Software Foundation; either version 2 # of the License, or (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 59 Temple Place - Suite 330, # Boston, MA 02111-1307, USA. # ################################################################## ###Last modification : 2/22/04 ###Preconditions of the main function "gbme" ##objects with no defaults #Y : a square matrix (n*n) ##objects with defaults #Xd : an n*n*rd array representing dyad-specific predictors #Xs : an n*rs matrix of sender-specific predictors #Xr : an n*rr matrix of receiver-specific predictors # fam : "gaussian" "binomial", or "poisson" # if fam="binomial" then an (n*n) matrix N of trials is needed # # pi.Sab = (s2a0,sab0,s2b0,v0) : parameters of inverse wishart dist # pi.s2u = vector of length 2: parameters of an inv-gamma distribution # pi.s2v = vector of length 2: parameters of an inv-gamma distribution # pi.s2z = k*2 matrix : parameters of k inv-gamma distributions # # pim.bd : vector of length rd, the prior mean for beta.d # piS.bd : a rd*rd pos def matrix, the prior covariance for beta.d # # pim.b0sr : vector of length 1+rs+rr, prior mean for (beta0,beta.s,beta.r) # piS.b0sr : a (1+rs+rr)*(1+rs+rr) pos def matrix, # the prior variance for (beta0,beta.s,beta.r) # k : dimension of latent space # directed : T or F. Is the data directed or undirected? If set to F, # estimation proceeds assuming Y[i,j]==Y[j,i], # and any supplied value of Xr is not used # also, only s2a0 and v0 in pi.Sab need be supplied, # and pi.s2v is not used. rr is zero, so pim.b0sr is # length 1+rs # starting values # beta.d #a vector of length rd # beta.u #a vector of length 1+rs+rr # s #a vector of length n # r #a vector of length n # z #an n*k matrix ###the main function gbme<-function( #data Y, Xd=array(dim=c(n,n,0)), Xs=matrix(nrow=dim(Y)[1],ncol=0), Xr=matrix(nrow=dim(Y)[1],ncol=0), #model specification fam="gaussian", N=NULL, k=0, directed=T, #compute starting values and emp Bayes for priors? estart=T, #priors - provide if estart=F pi.Sab=NULL, pi.s2u=NULL, pi.s2v=NULL, pi.s2z=NULL, pim.bd=NULL, piS.bd=NULL, pim.b0sr=NULL, piS.b0sr=NULL, #starting values - provide if estart=F beta.d=NULL, beta.u=NULL, s=NULL, r=NULL, z=NULL, #random seed, seed=0, #details of output NS=100000, #number of scans of mcmc to run odens=100, #output density owrite=T, #write output to a file? ofilename="OUT", #name of output file oround=5, #rounding of output zwrite=(k>0), #write z to file zfilename="Z", #outfile for z efilename="U", #outfile for z ffilename="V", #outfile for z sdigz=3, #rounding for z awrite=F,bwrite=F, afilename="A", bfilename="B" ) { ###set seed set.seed(seed) ###dimensions of everything n<<-dim(Y)[1] rd<<-dim(Xd)[3] rs<<-dim(Xs)[2] rr<<-dim(Xr)[2] diag(Y)<<-rep(0,n) ### ###column names for output file if(directed==T) { cnames<-c("k","scan","ll",paste("bd",seq(1,rd,length=rd),sep="")[0:rd],"b0", paste("bs",seq(1,rs,length=rs),sep="")[0:rs], paste("br",seq(1,rr,length=rr),sep="")[0:rr], "s2a","sab","s2b","s2e","rho", paste("s2e",seq(1,k,length=k),sep="")[0:k], paste("s2f",seq(1,k,length=k),sep="")[0:k] ) } if(directed==F) { cnames<-c("k","scan","ll",paste("bd",seq(1,rd,length=rd),sep="")[0:rd],"b0", paste("bs",seq(1,rs,length=rs),sep="")[0:rs], "s2a","s2e",paste("s2z",seq(1,k,length=k),sep="")[0:k] ) } ###design matrix for unit-specific predictors if(directed==T) { X.u<<-rbind( cbind( rep(.5,n), Xs, matrix(0,nrow=n,ncol=rr)), cbind( rep(.5,n), matrix(0,nrow=n,ncol=rs), Xr ) ) } if(directed==F) { X.u<<-cbind( rep(.5,n), Xs ) } ### ###construct an upper triangular matrix (useful for later computations) tmp<-matrix(1,n,n) tmp[lower.tri(tmp,diag=T)]<-NA UT<<-tmp ### ###get starting values and empirical bayes priors using a glm-type approach if(estart==T) {tmp<-gbme.glmstart(Y,Xd,Xs,Xr,fam,k,directed) priors<-tmp$priors ; startv<-tmp$startv } if(directed==T) { pi.Sab<-priors$pi.Sab; pi.s2u<-priors$pi.s2u; pi.s2v<-priors$pi.s2v pi.s2z<-priors$pi.s2z ; pim.bd<-priors$pim.bd ; piS.bd<-priors$piS.bd pim.b0sr<-priors$pim.b0sr ; piS.b0sr<-priors$piS.b0sr beta.d<-startv$beta.d ; beta.u<-startv$beta.u s<-startv$s ; r<-startv$r ; z<-startv$z ; e<-z ; f<-z } ### if(directed == F) { rho<-1 ;sv<-0 ; Xr<-Xs ; rr<-rs pi.s2u<-priors$pi.s2u ; pi.s2z<-priors$pi.s2z pim.bd<-priors$pim.bd ; piS.bd<-priors$piS.bd beta.d<-startv$beta.d ; beta.u<-startv$beta.u ; z<-startv$z s<-startv$s ; pi.s2a<-priors$pi.s2a ; pim.b0s<-priors$pim.b0s ; piS.b0s<-priors$piS.b0s Xr<-NULL ; rr<-0 ; r<-s*0 } ###Matrices for ANOVA on Xd, s, r tmp<-TuTv(n) #unit level effects design matrix for u=yij+yji, v=yij-yji Tu<-tmp$Tu Tv<-tmp$Tv if(directed==F) { Tu<-2*Tu[,1:n] } tmp<-XuXv(Xd) #regression design matrix for u=yij+yji, v=yij-yji Xu<-tmp$Xu Xv<-tmp$Xv XTu<<-cbind(Xu,Tu) XTv<<-cbind(Xv,Tv) tXTuXTu<<-t(XTu)%*%XTu tXTvXTv<<-t(XTv)%*%XTv ### ###redefining hyperparams if(directed==T){ Sab0<-matrix( c(pi.Sab[1],pi.Sab[2],pi.Sab[2],pi.Sab[3]),nrow=2,ncol=2) v0<-pi.Sab[4] } ### ###initializing error matrix E<-matrix(0,nrow=n,ncol=n) ### ###if using the linear link, then theta=Y if(fam=="gaussian"){ theta<-Y } ### ##### The Markov chain Monte Carlo algorithm for(ns in 1:NS){ ###update value of theta if(fam!="gaussian"){theta<- theta.betaX.d.srE.ef(beta.d,Xd,s,s*(1-directed)+r*directed,E,e,f) } ###impute any missing values if gaussian if(fam=="gaussian" & any(is.na(Y))) { rho<-(pi.s2u[2]-pi.s2v[2])/(pi.s2u[2]+pi.s2v[2]) se<-(pi.s2u[2]+pi.s2v[2])/4 mu<-theta.betaX.d.srE.ef(beta.d,Xd,s,s*(1-directed)+r*directed,E*0,e,f) imiss<-(1:n)[ is.na(Y%*%rep(1,n)) ] for(i in imiss){ for(j in (1:n)[is.na(Y[i,])]) { if( is.na(Y[i,j]) & is.na(Y[j,i]) ) { smp<-rmvnorm(c(mu[i,j],mu[j,i]),matrix( se*c(1,rho,rho,1),nrow=2,ncol=2)) theta[i,j]<-smp[1] ; theta[j,i]<-smp[2] } else{theta[i,j]<-rmvnorm( mu[i,j]+rho*(theta[j,i]-mu[j,i]), sqrt(se*(1-rho^2)))} } } } ###Update regression part tmp<-uv.E(theta-e%*%t(f)) #the regression part u<-tmp$u #u=yij+yji, i0){ s2e<-1/rgamma(k,pi.s2z[,1]+n/2,pi.s2z[,2]+diag( t(e)%*%e)/2) s2f<-1/rgamma(k,pi.s2z[,1]+n/2,pi.s2z[,2]+diag( t(f)%*%f)/2) #Gibbs for zs, using regression res<-theta-theta.betaX.d.srE.ef(beta.d,Xd,s,s*(1-directed)+r*directed,0*E, 0*e,0*f) s2u.res<-2*se*(1+rho) s2v.res<-2*se*(1-rho) for(i in sample(1:n)) { u.res<-(res[i,-i]+res[-i,i]) v.res<-(res[i,-i]-res[-i,i]) alp.j<-cbind( f[-i,],e[-i,] ) gam.j<-cbind( f[-i,],-e[-i,]) Sef<-chol2inv(chol( diag( 1/c(s2e,s2f),nrow=2*k) + t(alp.j)%*%alp.j/s2u.res + t(gam.j)%*%gam.j/s2v.res )) muef<-Sef%*%( t(alp.j)%*%u.res/s2u.res + t(gam.j)%*%v.res/s2v.res ) ef<-t(rmvnorm(muef,Sef)) e[i,]<-ef[1,1:k] ; f[i,]<-ef[1,k+1:k] } } ### if(k==0){sz<-NULL } #update E if( fam!="gaussian"){ E<-theta-theta.betaX.d.srE.ef(beta.d,Xd,s,s*(1-directed)+r*directed,E*0,e,f) #current E|beta.d,Xd,s,r,z UU<-VV<-Y*0 u<-rnorm(n*(n-1)/2,0,sqrt(su)) v<-rnorm(n*(n-1)/2,0,sqrt(sv)) UU[!is.na(UT)]<-u tmp<-t(UU) tmp[!is.na(UT)]<-u UU<-t(tmp) VV[!is.na(UT)]<-v tmp<-t(VV) tmp[!is.na(UT)]<--v VV<-t(tmp) Ep<-(UU+VV)/2 theta.p<-theta-E+Ep diag(theta.p)<-diag(theta)<-NA ##cases if(fam=="poisson") { mup<-exp(theta.p) ; mu<-exp(theta) M.lhr<-dpois(Y,mup,log=T) - dpois(Y,mu,log=T) } if(fam=="binomial"){ mup<-exp(theta.p)/(1+exp(theta.p)) mu<-exp(theta)/(1+exp(theta)) M.lhr<-dbinom(Y,N,mup,log=T) - dbinom(Y,N,mu,log=T) } M.lhr[ is.na(M.lhr) ]<-0 M.lhr<-(M.lhr+t(M.lhr))/( 2^(1-directed) ) RU<-matrix(log(runif(n*n)),nrow=n,ncol=n) RU[is.na(UT)]<-0 ; RU<-RU+t(RU) E[M.lhr>RU]<-Ep[M.lhr>RU] } ####output if(ns%%odens==0){ #compute loglikelihoods of interest if(fam=="poisson") { lpy.th<-sum( dpois(Y,exp(theta),log=T),na.rm=T) } if(fam=="binomial"){ lpy.th<-sum( dbinom(Y,N,exp(theta)/(1+exp(theta)), log=T),na.rm=T) } if(fam=="gaussian"){ E<-theta-theta.betaX.d.srE.ef(beta.d,Xd,s,s*(1-directed)+r*directed,E*0,e,f) tmp<-uv.E(theta-e%*%t(f)) #the regression part u<-tmp$u #u=yij+yji, iodens) { write.table(t(out),file=ofilename,append=T,quote=F, row.names=F,col.names=F) } } if(owrite==F) {cat(out,"\n") } if(zwrite==T) { write.table(signif(e,sdigz),file=efilename, append=T*(ns>odens),quote=F,row.names=F,col.names=F) } if(zwrite==T) { write.table(signif(f,sdigz),file=ffilename, append=T*(ns>odens),quote=F,row.names=F,col.names=F) } if(awrite==T){write.table(round(t(a),oround),file=afilename,append=T*(ns>odens),quote=F, row.names=F,col.names=F) } if(bwrite==T & directed==T){write.table(round(t(b),oround), file=bfilename,append=T*(ns>odens),quote=F, row.names=F,col.names=F) } } }} #####################End of main function : below are helper functions #### TuTv<-function(n){ Xu<-Xv<-NULL for(i in 1:(n-1)){ tmp<-tmp<-NULL if( i >1 ){ for(j in 1:(i-1)){ tmp<-cbind(tmp,rep(0,n-i)) } } tmp<-cbind(tmp,rep(1,n-i)) tmpu<-cbind(tmp,diag(1,n-i)) ; tmpv<-cbind(tmp,-diag(1,n-i)) Xu<-rbind(Xu,tmpu) ; Xv<-rbind(Xv,tmpv) } list(Tu=cbind(Xu,Xu),Tv=cbind(Xv,-Xv)) } #### XuXv<-function(X){ Xu<-Xv<-NULL if(dim(X)[3]>0){ for(r in 1:dim(X)[3]){ xu<-xv<-NULL for(i in 1:(n-1)){ for(j in (i+1):n){ xu<-c(xu,X[i,j,r]+X[j,i,r]) xv<-c(xv,X[i,j,r]-X[j,i,r]) }} Xu<-cbind(Xu,xu) Xv<-cbind(Xv,xv) } } list(Xu=Xu,Xv=Xv)} ### uv.E<-function(E){ u<- c( t( ( E + t(E) ) *UT ) ) u<-u[!is.na(u)] v<-c( t( ( E - t(E) ) *UT ) ) v<-v[!is.na(v)] list(u=u,v=v)} #### #### theta.betaX.d.srE.z<-function(beta.d,X.d,s,r,E,z){ m<-dim(X.d)[3] mu<-matrix(0,nrow=length(s),ncol=length(s)) if(m>0){for(l in 1:m){ mu<-mu+beta.d[l]*X.d[,,l] }} tmp<-mu+re(s,r,E,z) diag(tmp)<-0 tmp} #### theta.betaX.d.srE.ef<-function(beta.d,X.d,s,r,E,e,f){ m<-dim(X.d)[3] mu<-matrix(0,nrow=length(s),ncol=length(s)) if(m>0){for(l in 1:m){ mu<-mu+beta.d[l]*X.d[,,l] }} tmp<-mu+reef(s,r,E,e,f) diag(tmp)<-0 tmp} #### re<-function(a,b,E,z){ n<-length(a) matrix(a,nrow=n,ncol=n,byrow=F)+matrix(b,nrow=n,ncol=n,byrow=T)+E+z%*%t(z) } #### reef<-function(a,b,E,e,f){ n<-length(a) matrix(a,nrow=n,ncol=n,byrow=F)+matrix(b,nrow=n,ncol=n,byrow=T)+E+e%*%t(f) } #### #### rbeta.d.sr.gibbs<-function(u,v,su,sv,piS.bd,mu,Sab){ del<-Sab[1,1]*Sab[2,2]-Sab[1,2]^2 iSab<-rbind(cbind( diag(rep(1,n))*Sab[2,2]/del ,-diag(rep(1,n))*Sab[1,2]/del), cbind( -diag(rep(1,n))*Sab[1,2]/del,diag(rep(1,n))*Sab[1,1]/del) ) rd<-dim(as.matrix(piS.bd))[1] if(dim(piS.bd)[1]>0){ cov.beta.sr<-matrix(0,nrow=rd,ncol=2*n) iS<-rbind(cbind(solve(piS.bd),cov.beta.sr), cbind(t(cov.beta.sr),iSab)) } else{iS<-iSab} Sig<-chol2inv( chol(iS + tXTuXTu/su + tXTvXTv/sv ) ) #this may have a closed form expression #theta<-Sig%*%( (t(XTu)%*%u)/su + (t(XTv)%*%v)/sv + iS%*%mu) theta<-Sig%*%( t( (u%*%XTu)/su + (v%*%XTv)/sv ) + iS%*%mu) beta.sr<-rmvnorm(theta,Sig) list(beta.d=beta.sr[(rd>0):rd],s=beta.sr[rd+1:n],r=beta.sr[rd+n+1:n]) } #### #### rse.beta.d.gibbs<-function(g0,g1,x,XTx,s,r,beta.d){ n<-length(s) 1/rgamma(1, g0+choose(n,2)/2,g1+.5*sum( (x-XTx%*%c(beta.d,s,r))^2 ) ) } #### #### rbeta.sr.gibbs<-function(s,r,X.u,pim.b0sr,piS.b0sr,Sab) { del<-Sab[1,1]*Sab[2,2]-Sab[1,2]^2 iSab<-rbind(cbind( diag(rep(1,n))*Sab[2,2]/del ,-diag(rep(1,n))*Sab[1,2]/del), cbind( -diag(rep(1,n))*Sab[1,2]/del,diag(rep(1,n))*Sab[1,1]/del) ) S<-solve( solve(piS.b0sr) + t(X.u)%*%iSab%*%X.u ) mu<-S%*%( solve(piS.b0sr)%*%pim.b0sr+ t(X.u)%*%iSab%*%c(s,r)) rmvnorm( mu,S) } #### #### rSab.gibbs<-function(a,b,S0,v0){ n<-length(a) ab<-cbind(a,b) Sn<-S0+ (t(ab)%*%ab) solve(rwish(solve(Sn),v0+n) ) } #### rmvnorm<-function(mu,Sig2){ R<-t(chol(Sig2)) R%*%(rnorm(length(mu),0,1)) +mu } #### rwish<-function(S0,nu){ S<-S0*0 for(i in 1:nu){ z<-rmvnorm(rep(0,dim(as.matrix(S0))[1]), S0) S<-S+z%*%t(z) } S } ### Procrustes transformation: rotation and reflection proc.rr<-function(Y,X){ k<-dim(X)[2] A<-t(Y)%*%( X%*%t(X) )%*%Y eA<-eigen(A,symmetric=T) Ahalf<-eA$vec[,1:k]%*%diag(sqrt(eA$val[1:k]),nrow=k)%*%t(eA$vec[,1:k]) t(t(X)%*%Y%*%solve(Ahalf)%*%t(Y)) } ### gbme.glmstart<-function(Y,Xd,Xs,Xr,fam,k,directed){ #generate starting values and emp bayes priors from #approximate mles y<-x<-NULL for (i in 1:n) { for (j in (1:n)[-i]) { y<-c(y,Y[i,j]) x<-rbind(x,c(i,j,Xd[i,j,])) }} if(rd>0) { fit1<-glm(y~-1+as.factor(x[,1])+as.factor(x[,2])+x[,-(1:2)],family=fam) } if(rd==0){ fit1<-glm(y~-1+as.factor(x[,1])+as.factor(x[,2]),family=fam) } s<-fit1$coef[1:n] r<-c(0,fit1$coef[n+1:(n-1)]) if(rd>0){ beta.d<-fit1$coef[(2*n):length(fit1$coef)] #beta.d } la<-lm( s~-1+cbind(rep(.5,n),Xs)); a<-la$res lb<-lm( r~-1+cbind(rep(.5,n),Xr)); b<-lb$res b0<-(la$coef[1]+lb$coef[1])/2 s1<-s+( b0-la$coef[1])/2 #s r1<-r+( b0-lb$coef[1])/2 #r la<-lm( s1~-1+cbind(rep(.5,n),Xs)); a<-la$res lb<-lm( r1~-1+cbind(rep(.5,n),Xr)); b<-lb$res beta.u<-c(lb$coef[1], la$coef[-1],lb$coef[-1]) Sab<-var(cbind(a,b)) #Sab Res<-matrix(0,nrow=n,ncol=n) l<-1 for( m in 1:dim(x)[1] ) { if( !is.na(y[m]) ) { Res[x[m,1],x[m,2]]<-fit1$res[l] ; l<-l+1 } } ###### Res[Res>quantile(Res,.95,na.rm=T) ] <- quantile(Res,.95,na.rm=T) diag(Res)<-rep(0,n) if(k>0){ for(i in 1:50){ tmp<-eigen( .5*(Res+t(Res) )) z<-tmp$vec[,1:k] %*%diag(sqrt(tmp$val[1:k]),nrow=k) diag(Res)<-diag(z%*%t(z)) } sz<-var(c(z)) #z,sz E<-Res-z%*%t(z) #E } if(k==0){E<-Res ; z<-matrix(0,nrow=n,ncol=1) ;sz<-.1} u<-E+t(E) u<-c(u[!is.na(UT)]) #su su<-var(u) #sv v<-E-t(E) v<-c(v[!is.na(UT)]) sv<-var(v) rho<-(su-sv)/(su+sv) se<-(su+sv)/4 ###construct priors centered on these empirical values pi.Sab<-c(Sab[1,1],Sab[2,1],Sab[2,2],4) pi.s2u<-c(2,su) pi.s2v<-c(2,sv) pi.s2z<-cbind( rep(2,k), diag( t(z)%*%z)/n ) if(rd>0){ pim.bd<-c(beta.d) piS.bd<-as.matrix((n*(n-1))^2*summary(fit1)$cov.unscaled[(2*n):length(fit1$coef), (2*n):length(fit1$coef)] ) } if(rd==0){ beta.d<-pim.bd<-NULL ; piS.bd<-matrix(nrow=0,ncol=0) } if(directed==T) { pim.b0sr<-beta.u piS.b0sr<-rbind(cbind( summary(la)$cov[-1,-1], matrix(0,nrow=rs,ncol=rr)), # cbind(matrix(0,nrow=rs,ncol=rr),summary(lb)$cov[-1,-1])) # error found 7-7-6 cbind(matrix(0,nrow=rr,ncol=rs),summary(lb)$cov[-1,-1])) piS.b0sr<-rbind(c(summary(la)$cov[1,-1],summary(lb)$cov[1,-1]),piS.b0sr) piS.b0sr<-cbind(c(summary(la)$cov[1,1:(rs+1)], summary(lb)$cov[1,-1]),piS.b0sr) piS.b0sr<-as.matrix(piS.b0sr*n*n) pim.b0s<-NULL ; piS.b0s<-NULL ;pi.s2a<-NULL } if(directed==F) { pim.b0sr<-NULL ; piS.b0sr<-NULL ; pi.Sab<-NULL s<-(s+r)/2 ; r<-s*0 tmp2<-lm(s~-1+cbind( rep(.5,n),Xs)) beta.u<-tmp2$coef ; pi.s2a<-c(2,var(tmp2$res)) ; pim.b0s<-tmp2$coef ; piS.b0s<-summary(tmp2)$cov*n*n } ### list( priors=list(pi.Sab=pi.Sab,pi.s2u=pi.s2u,pi.s2v=pi.s2v,pi.s2z=pi.s2z, pim.bd=pim.bd,piS.bd=piS.bd, pim.b0sr=pim.b0sr,piS.b0sr=piS.b0sr, pim.b0s=pim.b0s, piS.b0s=piS.b0s , pi.s2a=pi.s2a) , startv=list(beta.d=beta.d,beta.u=beta.u,s=s,r=r,z=z) ) } #### rbeta.d.s.gibbs<-function(u,su,piS.bd,mu,s2a){ iSa<-diag(rep(1/s2a,n)) if(dim(piS.bd)[1]>0){ cov.beta.s<-matrix(0,nrow=rd,ncol=n) iS<-rbind(cbind(solve(piS.bd),cov.beta.s), cbind(t(cov.beta.s),iSa)) } else{iS<-iSa} Sig<-chol2inv( chol(iS + tXTuXTu/su ) ) #this may have a closed form expression theta<-Sig%*%( t( (u%*%XTu)/su) + iS%*%mu) beta.s<-rmvnorm(theta,Sig) list(beta.d=beta.s[(rd>0):rd],s=beta.s[rd+1:n]) } #### #### rbeta.s.gibbs<-function(s,X.u,pim.b0s,piS.b0s,s2a) { iSa<-diag(rep(1/s2a,n)) S<-solve( solve(piS.b0s) + t(X.u)%*%iSa%*%X.u ) mu<-S%*%( solve(piS.b0s)%*%pim.b0s+ t(X.u)%*%iSa%*%s) rmvnorm( mu,S) }