R/bnpLongMulti.R

Defines functions plot3.bcmg postSigma2 postSigma lmrm

Documented in lmrm

# lmrm, postSigma, postSigma2, plot3.bcmg  
lmrm <- function(formula,
                 data=list(),centre=TRUE,
                 id,
                 time,
                 sweeps,burn=0,thin=1,seed,StorageDir,
                 c.betaPrior="IG(0.5,0.5*n*p)", 
                 pi.muPrior="Beta(1,1)", 
                 c.alphaPrior="IG(1.1,1.1)",
                 pi.phiPrior="Beta(1,1)", 
                 c.psiPrior="HN(2)",
                 sigmaPrior="HN(2)", 
                 pi.sigmaPrior="Beta(1,1)", 
                 corr.Model=c("common",nClust=1), 
                 DP.concPrior="Gamma(5,2)",
                 c.etaPrior="IG(0.5,0.5*samp)", 
                 pi.nuPrior="Beta(1,1)", 
                 pi.fiPrior="Beta(1,1)",
                 c.omegaPrior="IG(1.1,1.1)",
                 sigmaCorPrior="HN(2)",
                 tuneCalpha,
                 tuneSigma2,
                 tuneCbeta,
                 tuneAlpha,
                 tuneSigma2R,
                 tuneR,
                 tuneCpsi,
                 tuneCbCor,
                 tuneOmega,
                 tuneComega,
                 tau,FT=1,...){
    #Samples etc
    if (thin <= 0) thin <- 1
    thin <- as.integer(thin)
    sweeps <- as.integer(sweeps)
    burn <- as.integer(burn)
    if (missing(sweeps)) stop("provide sweeps argument")
    nSamples<-0
    LASTsw<-0
    if (sweeps > 0 && (sweeps-burn) > 0){
		nSamples <- length(seq(1,(sweeps-burn),by=thin))
		LASTsw<-seq(burn,by=thin,length.out=nSamples)[nSamples]
	}
	LASTWB<-1	
    # Match call
    call <- match.call(expand.dots = FALSE)
    call2 <- match.call.defaults()
    #Data
    if (!is.list(data) && !is.data.frame(data))
        data <- as.data.frame(data)
    #Formula & data dimensions
    p <- length(as.Formula(formula))[1]    
    if (length(as.Formula(formula))[2] != 5) stop("ambiguous definition of regression model")    
    formula.m<-formula(as.Formula(formula),lhs=0,rhs=1)
    formula.v<-formula(as.Formula(formula),lhs=0,rhs=2)
    formula.d<-formula(as.Formula(formula),lhs=0,rhs=3)
    formula.corm<-formula(as.Formula(formula),lhs=0,rhs=4)
    formula.corv<-formula(as.Formula(formula),lhs=0,rhs=5)
    # Responses, design matrices, indicators
    Y<-NULL
    varsY<-list()
    for (i in 1:p){
        trms<-terms.formula(formula(as.Formula(formula,~1),lhs=i,rhs=0))
        Y<-cbind(Y,with(data,eval(attr(trms,"variables")[[2]])))
        varsY[[i]] <- as.character(attr(trms,"variables")[[2]])
    }
    # Null Deviance
    nullDeviance <- 0
    for (i in 1:p) nullDeviance <- nullDeviance - 2 * logLik(lm(Y[,i] ~ 1))
    # Response etc
    if (missing(id)) stop("id needed")    
    id<-eval(substitute(id), data)
    if (is.null(id)) stop("unrecognised id")
    n<-length(unique(id))
    #print("n");print(n)
    niVec<-table(id)[1:n] 
    niMax<-max(niVec)
    N<-sum(niVec)
    cusumniVec<-c(0,cumsum(niVec))
    cusumC<-c(0,cumsum(niVec*(niVec-1)/2))     
    # Times
    varTime<-as.character(substitute(list(time))[-1])
    time2<-eval(substitute(time), data)
    LUT<-length(unique(time2)) 
    SUT<-sort(unique(time2))
    FUT<-table(time2)
    intime<-time2
    for (i in 1:length(intime)) intime[i] <- which(intime[i]==SUT)-1   
    intime2 <- rep(intime,each=p)
    #Desing matrix C (dependence)
    dataNew<-as.data.frame(cbind(data,lag=rnorm(dim(data)[1])))
    bb<-DM(formula=formula.d,data=dataNew,n=dim(data)[1],centre=centre)
    vars <- unlist(bb$vars)
    if ("lag" %in% vars)
        vars <- vars[-which(vars=="lag")]   
    C<-NULL
    for (i in 1:n) 
        if (niVec[i] > 1) C<-rbind(C,cbind(cusumniVec[i]+rep(seq(2,niVec[i],1),seq(1,niVec[i]-1,1)), 
                                           cusumniVec[i]+unlist(sapply(1:(niVec[i]-1), function(i) seq(1,i,1)))                                     
                                     ))
    lag<-time2[C[,1]]-time2[C[,2]]
    lag<-data.frame(lag=lag)    
    if (length(vars)>0){
        for (k in 1:length(vars)){
            V<-data[,vars[k]]
            lag<-data.frame(lag,V[C[,1]])
            colnames(lag)[k+1] <- vars[[k]]
	    }
	}
	#print(lag)
	if (niMax>1){    
        XYK<-DM(formula=formula.d,data=lag,n=dim(lag)[1],centre=centre)
        C<-as.matrix(XYK$X)
    }else{
    	XYK<-NULL
	    C<-NULL
	}
    Cknots<-XYK$Rknots
    storeMeanVectorC<-XYK$meanVector
    storeIndicatorC<-XYK$indicator
    LK<-NCOL(C)
    if (!is.null(C)) C<-t(C)
    vecLK<-table(XYK$assign)
    NK<-length(vecLK)
    cusumVecLK<-c(0,cumsum(vecLK))
    assignC<-XYK$assign
    labelsC<-XYK$labels
    countC<-XYK$count
    varsC<-XYK$vars
    is.Dc<-XYK$is.D
    which.SpecC<-XYK$which.Spec
    formula.termsC<-XYK$formula.terms
    #Desing matrix X (mean)
    XYK<-DM(formula=formula.m,data=data,n=N,centre=centre)
    X<-as.matrix(XYK$X)
    Xknots<-XYK$Rknots
    storeMeanVectorX<-XYK$meanVector
    storeIndicatorX<-XYK$indicator
    LG<-NCOL(X)-1
    vecLG<-table(XYK$assign)[-1]
    NG<-length(vecLG)
    cusumVecLG<-c(0,cumsum(vecLG))
    assignX<-XYK$assign
    labelsX<-XYK$labels
    countX<-XYK$count
    varsX<-XYK$vars
    is.Dx<-XYK$is.D
    which.SpecX<-XYK$which.Spec
    formula.termsX<-XYK$formula.terms    
    #print("2")
    #Desing matrix Z (variances)
    ZK<-DM(formula=formula.v,data=data,n=N,centre=centre)
    Z<-as.matrix(ZK$X)
    Zknots<-ZK$Rknots
    storeMeanVectorZ<-ZK$meanVector
    storeIndicatorZ<-ZK$indicator
    LD<-NCOL(Z)-1
    labelsZ<-ZK$labels
    countZ<-ZK$count
    varsZ<-ZK$vars
    which.SpecZ<-ZK$which.Spec
    formula.termsZ<-ZK$formula.terms
    is.Dz<-ZK$is.D
    assignZ<-ZK$assign    
    #vecLD<-table(assignZ)[-1]
    #cusumVecLD<-c(0,cumsum(vecLD))
    ND<-max(assignZ)#length(vecLD)    
    repeats<-rep(1,ND)
    repeats[which(is.Dz==1)]<-table(assignZ)[which(is.Dz==1)+1]
    oneHE<-1; if (ND > 0) oneHE<-rep(c(1:ND),repeats)
    assignZ2<-ZK$assign2        
    vecLD<-table(assignZ2)[-1]
    cusumVecLD<-c(0,cumsum(vecLD))
    ND2<-max(assignZ2)#length(vecLD)    
    isDz2<-as.numeric(is.Dz)[oneHE]
    if (is.na(isDz2)) isDz2<-1
    MVLD <- 1
    if (LD > 0 || LK > 0) MVLD<-max(vecLD,vecLK)
    #Design matrix Xc (mean of correlations)
    d<-p*(p-1)/2
    t <- rep(SUT,each=d)
    dataCor<-data.frame(t)
    if (p > 1){
		lab<-DM(formula=formula.corm,data=data,n=1,centre=centre)        
        if (length(lab$vars) > 0) colnames(dataCor)<-lab$vars[[1]]   
		XYK<-DM(formula=formula.corm,data=dataCor,n=(LUT*d),centre=centre)
		Xc<-as.matrix(XYK$X)
	}else{
	    XYK<-NULL
	    Xc<-NULL
	}    
    Xcknots<-XYK$Rknots
    storeMeanVectorXc<-XYK$meanVector
    storeIndicatorXc<-XYK$indicator
    LGc<-NCOL(Xc)-1
    vecLGc<-table(XYK$assign)[-1]
    NGc<-length(vecLGc)
    cusumVecLGc<-c(0,cumsum(vecLGc))
    assignXc<-XYK$assign
    labelsXc<-XYK$labels
    countXc<-XYK$count
    varsXc<-XYK$vars
    is.Dxc<-XYK$is.D
    which.SpecXc<-XYK$which.Spec
    formula.termsXc<-XYK$formula.terms
    #Design matrix Zc (variance of correlations)
    if (p > 1){
		XYK<-DM(formula=formula.corv,data=dataCor,n=(LUT*d),centre=centre)    
        Zc<-as.matrix(XYK$X)
    }else{
		XYK<-NULL
		Zc<-NULL
	}
    Zcknots<-XYK$Rknots
    storeMeanVectorZc<-XYK$meanVector
    storeIndicatorZc<-XYK$indicator
    LDc<-NCOL(Zc)-1
    vecLDc<-table(XYK$assign)[-1]
    MVLD<-max(vecLDc,MVLD)
    NDc<-length(vecLDc)
    cusumVecLDc<-c(0,cumsum(vecLDc))
    assignZc<-XYK$assign
    labelsZc<-XYK$labels
    countZc<-XYK$count
    varsZc<-XYK$vars
    is.Dzc<-XYK$is.D
    which.SpecZc<-XYK$which.Spec
    formula.termsZc<-XYK$formula.terms   
    #Initialize covariance & correlation matrix
    #print("3")
    LASTR<-LASTD<-LASTE<-rep(1,LUT)    
    if (p>1){        
        LASTR<-LASTD<-LASTE<-NULL
        for (t in SUT){
            Res<-NULL
            for (i in 1:p){            
				#print(head(X))
				#print(grep("(",colnames(X),fixed=TRUE)) 
				#RM<-grep("(",colnames(X),fixed=TRUE)
				#print(RM)    			
                RM<-grep("(",colnames(X),fixed=TRUE)[-1]
				#print(RM)
				Xinit<-X[time2==t,]
                if (length(RM) > 0) Xinit<-Xinit[,-RM]
                #print(head(Xinit))
                lm1<-lm(Y[time2==t,i] ~ Xinit)
                Res<-cbind(Res,residuals(lm1))
			}
            CR<-0.9*cov(Res)+0.1*diag(mean(eigen(cov(Res))$values),p)         
            LASTR<-c(LASTR,c(cov2cor(CR)))            
            D<-matrix(0,p,p)
            diag(D)<-sqrt(diag(CR))
            LASTD<-c(LASTD,c(D))
            LASTE<-c(LASTE,c(CR))
            #LASTmuR<-mean(LASTR[upper.tri(LASTR)])
            #LASTsigma2R<-1
            #if (p > 2) LASTsigma2R<-var(LASTR[upper.tri(LASTR)])
            #print(CR)
            #print(D)
            #print(cov2cor(CR))            
		}
    }  
    LASTAll<-c(LASTR,LASTD,LASTE)    
    #print("4")
    #print(LASTAll)
    #print("2")     
    #Prior for pi.mu
    if (!length(pi.muPrior)==1 && !length(pi.muPrior)==NG && !length(pi.muPrior)==(p*NG))
        stop("pi.muPrior has incorrect dimension")
    pimu<-NULL
    for (k in 1:length(pi.muPrior)){
        sp<-strsplit(pi.muPrior[k],"Beta\\(")
        sp<-strsplit(sp[[1]][2],"\\)")
        sp<-strsplit(sp[[1]][1],",")
        pimu<-c(pimu,as.numeric(sp[[1]]))
    }
    if (length(pi.muPrior)==1) pimu<-rep(pimu,p*NG)
    if (length(pi.muPrior)==NG) pimu<-rep(pimu,p)
    #Prior for pi.phi
    if (!length(pi.phiPrior)==1 && !length(pi.phiPrior)==NK && !length(pi.phiPrior)==(p*p*NK))
        stop("pi.phiPrior has incorrect dimension")
    piphi<-NULL
    for (k in 1:length(pi.phiPrior)){
        sp<-strsplit(pi.phiPrior[k],"Beta\\(")
        sp<-strsplit(sp[[1]][2],"\\)")
        sp<-strsplit(sp[[1]][1],",")
        piphi<-c(piphi,as.numeric(sp[[1]]))
    }
    if (length(pi.phiPrior)==1) piphi<-rep(piphi,p*p*NK)
    if (length(pi.phiPrior)==NK) piphi<-rep(piphi,p*p)
    #Prior for c.beta
    if (!length(c.betaPrior)==1)
        stop("c.betaPrior has incorrect dimension")
    #c.betaPrior<-sub("samp","p*n",c.betaPrior)
    sp<-strsplit(c.betaPrior,"IG\\(")
    sp<-strsplit(sp[[1]][2],"\\)")
    sp<-strsplit(sp[[1]][1],",")
    cetaParams<-c(as.numeric(sp[[1]][1]),eval(parse(text=sp[[1]][2])))
    #print(cetaParams)
    #Prior for c.eta
    if (!length(c.etaPrior)==1)
        stop("c.etaPrior has incorrect dimension")
    c.etaPrior<-sub("samp","LUT*d",c.etaPrior)
    sp<-strsplit(c.etaPrior,"IG\\(")
    sp<-strsplit(sp[[1]][2],"\\)")
    sp<-strsplit(sp[[1]][1],",")
    cetaCorParams<-c(as.numeric(sp[[1]][1]),eval(parse(text=sp[[1]][2])))
    #Prior for pi.sigma
    if (ND > 0){
        if (!length(pi.sigmaPrior)==1 && !length(pi.sigmaPrior)==ND && !length(pi.sigmaPrior)==(p*ND))
            stop("pi.sigmaPrior has incorrect dimension")
        if (length(pi.sigmaPrior)==(p*ND)) 
            pi.sigmaPrior<-pi.sigmaPrior[oneHE+rep(seq(from=0,to=p*ND-1,by=ND),each=length(oneHE))]
        if (length(pi.sigmaPrior)==ND) pi.sigmaPrior<-pi.sigmaPrior[oneHE]    
        pisigma<-NULL
        for (k in 1:length(pi.sigmaPrior)){
            sp<-strsplit(pi.sigmaPrior[k],"Beta\\(")
            sp<-strsplit(sp[[1]][2],"\\)")
            sp<-strsplit(sp[[1]][1],",")
            pisigma<-c(pisigma,as.numeric(sp[[1]]))
        }
        if (length(pi.sigmaPrior)==1) pisigma<-rep(pisigma,p*ND2)
        if (length(pi.sigmaPrior)==ND2) pisigma<-rep(pisigma,p)
	}else{pisigma<-1}
    #Prior for c.alpha
    if (!length(c.alphaPrior)==1 && !length(c.alphaPrior)==p)
        stop("c.alphaPrior has incorrect dimension")
    specials<-c("HN","IG")
    calphaParams<-NULL
    HNca<-vector()
    for (k in 1:length(c.alphaPrior)){
        sp<-strsplit(c.alphaPrior[k],"\\(")
        if (sp[[1]][1] %in% specials){
            if (match(sp[[1]][1],specials)==1) HNca[k]<-1
            if (match(sp[[1]][1],specials)==2) HNca[k]<-0
        } else stop("unrecognised prior for c.alpha")
        sp<-strsplit(sp[[1]][2],"\\)")
        sp<-strsplit(sp[[1]][1],",")
        calphaParams<-c(calphaParams,as.numeric(sp[[1]]))
    }
    if (length(c.alphaPrior)==1){
        calphaParams<-rep(calphaParams,p)
        HNca<-rep(HNca,p)
	}
	#Prior for cPsi
    if (!length(c.psiPrior)==1 && !length(c.psiPrior)==(p*p))
        stop("c.psiPrior has incorrect dimension")
    specials<-c("HN","IG")
    cpsiParams<-NULL
    HNcpsi<-vector()
    for (k in 1:length(c.psiPrior)){
        sp<-strsplit(c.psiPrior[k],"\\(")
        if (sp[[1]][1] %in% specials){
            if (match(sp[[1]][1],specials)==1) HNcpsi[k]<-1
            if (match(sp[[1]][1],specials)==2) HNcpsi[k]<-0
        } else stop("unrecognised prior for c.psi")
        sp<-strsplit(sp[[1]][2],"\\)")
        sp<-strsplit(sp[[1]][1],",")
        cpsiParams<-c(cpsiParams,as.numeric(sp[[1]]))
	}
	if (length(c.psiPrior)==1){
        cpsiParams<-rep(cpsiParams,(p*p))
        HNcpsi<-rep(HNcpsi,(p*p))
	}    
    #Prior for sigma2_{zero k}
    if (!length(sigmaPrior)==1 && !length(sigmaPrior)==p)
        stop("sigmaPrior has incorrect dimension")
    specials<-c("HN","IG")
    sigmaParams<-NULL
    HNsg<-vector()
    for (k in 1:length(sigmaPrior)){
        sp<-strsplit(sigmaPrior[k],"\\(")
        if (sp[[1]][1] %in% specials){
            if (match(sp[[1]][1],specials)==1) HNsg[k]<-1
            if (match(sp[[1]][1],specials)==2) HNsg[k]<-0
        } else stop("unrecognised prior for sigma^2")
        sp<-strsplit(sp[[1]][2],"\\)")
        sp<-strsplit(sp[[1]][1],",")
        sigmaParams<-c(sigmaParams,as.numeric(sp[[1]]))
    }
    if (length(sigmaPrior)==1){
        sigmaParams<-rep(sigmaParams,p)
        HNsg<-rep(HNsg,p)
	}
    #Prior for pi.nu
    if (!length(pi.nuPrior)==NGc) if (!length(pi.nuPrior)==1) 
        stop("pi.nuPrior has incorrect dimension")
    pinu<-NULL
    for (k in 1:length(pi.nuPrior)){
        sp<-strsplit(pi.nuPrior[k],"Beta\\(")
        sp<-strsplit(sp[[1]][2],"\\)")
        sp<-strsplit(sp[[1]][1],",")
        pinu<-c(pinu,as.numeric(sp[[1]]))
    }
    if (length(pi.nuPrior)==1) pinu<-rep(pinu,NGc)
    #Prior for pi.fi
    if (!length(pi.fiPrior)==NDc) if (!length(pi.fiPrior)==1) 
        stop("pi.fiPrior has incorrect dimension")
    pifi<-NULL
    for (k in 1:length(pi.fiPrior)){
        sp<-strsplit(pi.fiPrior[k],"Beta\\(")
        sp<-strsplit(sp[[1]][2],"\\)")
        sp<-strsplit(sp[[1]][1],",")
        pifi<-c(pifi,as.numeric(sp[[1]]))
    }
    if (length(pi.fiPrior)==1) pifi<-rep(pifi,ND)    
    #Prior for c.omega
    if (!length(c.omegaPrior)==1)
        stop("c.omegaPrior has incorrect dimension")
    specials<-c("HN","IG")
    sp<-strsplit(c.omegaPrior,"\\(")
    if (sp[[1]][1] %in% specials){ 
        if (match(sp[[1]][1],specials)==1) HNco<-1
        if (match(sp[[1]][1],specials)==2) HNco<-0
    } else stop("unrecognised prior for c.omega")
    sp<-strsplit(sp[[1]][2],"\\)")
    sp<-strsplit(sp[[1]][1],",")
    comegaParams<-as.numeric(sp[[1]])
	#Prior for sigma2Cor
	if (!length(sigmaCorPrior)==1)
        stop("sigmaCorPrior has incorrect dimension")
    HNscor<-0    
    specials<-c("HN","IG")
    sp<-strsplit(sigmaCorPrior,"\\(")
    if (sp[[1]][1] %in% specials){ 
        if (match(sp[[1]][1],specials)==1 && p > 1) HNscor<-1
        if (match(sp[[1]][1],specials)==2) HNscor<-0
    } else stop("unrecognised prior for sigma2Cor")
    sp<-strsplit(sp[[1]][2],"\\)")
    sp<-strsplit(sp[[1]][1],",")
    sigmaCorParams<-as.numeric(sp[[1]])
    #Cor model
    corModels<-c("common","groupC","groupV")
    mcm<-4+match(corr.Model[1],corModels)
    if (p==1 || p==2) mcm <- 5
    if (is.na(mcm)) stop("unrecognised correlation model")
    H <- G <- 1
    if (mcm==6){
        H<-as.numeric(corr.Model[2])
        if (H==1) {mcm<-5; warning("Common correlations model specified with nClust = 1")}
        if (is.na(H) || (!H%%1==0) || H==0) {H <- p*(p-1)/2; warning(cat("mispecified number of clusters. nClust set to ",H,"\n"))}
    }
    if (mcm==7){
        G<-as.numeric(corr.Model[2])
        if (G==1) {mcm<-5; warning("Common correlations model specified with nClust = 1")}
        if (is.na(G) || (!G%%1==0) || G==0) {G <- p; warning(cat("mispecified number of clusters. nClust set to ",G,"\n"))}
        H<-G*(G-1)/2+G #min(d,G*(G-1)/2+G) #min(G,abs(p-G))
	}
    #Prior for alpha DP
    if (mcm == 6 || mcm==7){
        sp<-strsplit(DP.concPrior,"Gamma\\(")
        sp<-strsplit(sp[[1]][2],"\\)")
        sp<-strsplit(sp[[1]][1],",")
        DPparams<-as.numeric(sp[[1]])
        DPparams <- c(DPparams,0.01)
    }
    #Seed
    if (missing(seed)) seed<-as.integer((as.double(Sys.time())*1000+Sys.getpid()) %% 2^31)
    # Storage directory & files
    if (!missing(StorageDir)){
        StorageDir <- path.expand(StorageDir)
        ncharwd <- nchar(StorageDir)}
    if (!missing(StorageDir)) if (!(substr(StorageDir,ncharwd,ncharwd)=="/")) StorageDir <- paste(StorageDir,"/",sep="")
    if (!missing(StorageDir)) if (!dir.exists(StorageDir)) dir.create(StorageDir, recursive = TRUE)
    if (missing(StorageDir)) stop("provide a storage directory via argument StorageDir")
    FL <- c("gamma", "cbeta", "delta", "alpha", "sigma2", "calpha", "beta", "psi", "ksi", "cpsi", "R", "nu", "fi", "omega",
            "sigma2R","ceta","comega","eta", "deviance", 
            "compAlloc", "nmembers", "DPconc",
            "compAllocV","nmembersV",
            "DE", "test","nu.ls","eta.ls","nmembers.ls","clusters","probs","tune",
            "muR", "phi2")
    for (i in 1:length(FL)){
        oneFile <- paste(StorageDir, paste("BNSP",FL[i], "txt",sep="."),sep="/")
        if (file.exists(oneFile)) file.remove(oneFile)
	}
    #Tuning Parameters            
    if (missing(tuneCalpha)) tuneCalpha<-rep(1,p)
    if (!length(tuneCalpha)==p) tuneCalpha<-rep(mean(tuneCalpha),p)        
    if (missing(tuneSigma2)) tuneSigma2<-rep(1,p)
    if (!length(tuneSigma2)==p) tuneSigma2<-rep(mean(tuneSigma2),p)                 
    if (missing(tuneCbeta)) tuneCbeta<-100
    if (ND>0){
        if (missing(tuneAlpha) || !length(tuneAlpha)==ND || !length(tuneAlpha)==(ND*p)){
		    tuneAlpha<-rep(5,ND)
		    tuneAlpha[which(is.Dz==1)]<-1
	    }    
	    if (length(tuneAlpha)==(ND*p)) 
            tuneAlpha<-tuneAlpha[oneHE+rep(seq(from=0,to=p*ND-1,by=ND),each=length(oneHE))]                
        if (length(tuneAlpha)==ND) tuneAlpha<-tuneAlpha[oneHE]
        if (length(tuneAlpha)==ND2) tuneAlpha<-rep(tuneAlpha,p)
    }else{tuneAlpha<-1}
    
    if (missing(tuneSigma2R)) tuneSigma2R<-1
    if (missing(tuneR)) tuneR<-rep(40*(p+2)^3,LUT)
    tuneR[which(tuneR<p+2)]<-p+2
    if (!length(tuneR)==LUT) tuneR<-rep(mean(tuneR),LUT)
    if (missing(tuneCpsi)) tuneCpsi<-rep(5,p*p)
    if (!length(tuneCpsi)==(p*p)) tuneCpsi<-rep(mean(tuneCpsi),p*p)
	if (missing(tuneCbCor)) tuneCbCor<-10
	if (missing(tuneOmega)) tuneOmega<-rep(5,NDc)
	if (!length(tuneOmega)==NDc) tuneOmega<-rep(mean(tuneOmega),NDc)
    if (missing(tuneComega)) tuneComega<-1
	if (missing(tau)) tau = 0.01
    #Block size selection
    #if (missing(blockSizeProbG)){
    blockSizeProbG <- rep(0,LG)
    blockSizeProbG[1:5]<-c(10,25,30,25,10)
	#}
    #if (missing(blockSizeProbD)){
	blockSizeProbD <- rep(0,LD)
    blockSizeProbD[1:5]<-c(10,25,30,25,10)
    #}
    #if (missing(blockSizeProbK)){
	blockSizeProbK <- rep(0,LK)
    blockSizeProbK[1:5]<-c(10,25,30,25,10)
    #}
    #if (missing(blockSizeProbGc)){ 
    blockSizeProbGc <- rep(0,LGc) 
    blockSizeProbGc[1:5]<-c(10,25,30,25,10)
	#}
    #if (missing(blockSizeProbD)){
	blockSizeProbDc <- rep(0,LDc)
    blockSizeProbDc[1:5]<-c(10,25,30,25,10)
    #}
    maxBSG <- max(which(blockSizeProbG>0))
    maxBSD <- max(which(blockSizeProbD>0))
    maxBSK <- max(which(blockSizeProbK>0))
    maxBSGc <- max(which(blockSizeProbGc>0))
    maxBSDc <- max(which(blockSizeProbDc>0))
    
    #print(as.double(c(blockSizeProbG, blockSizeProbD, blockSizeProbK, blockSizeProbGc, blockSizeProbDc)))
    #print(as.integer(c(maxBSG,maxBSD,maxBSK,maxBSGc,maxBSDc)))
    
    #Deviance
    deviance <- c(0,0)
    #Cont
    cont <- 0  
    #tol
    tol <- sqrt(.Machine$double.eps)  
    #Call C
    if (mcm==5){
        out<-.C("longmult",
        as.integer(seed), as.character(StorageDir), 
        as.integer(c(sweeps,burn,thin,n,p,N,MVLD,mcm)),
        as.integer(niVec), as.integer(cusumniVec), as.integer(intime2), as.integer(intime), 
        as.integer(c(niMax,LUT)), as.integer(FUT),
        as.integer(cusumC), as.double(C), as.double(c(t(Y))), as.double(t(X)), as.double(Z), 
        as.double(Xc),as.double(Zc),
        as.integer(c(LG,LD,LK,LGc,LDc)), as.integer(c(NG,ND2,NK,NGc,NDc)), 
        as.integer(vecLG), as.integer(vecLD), as.integer(vecLK), as.integer(vecLGc),as.integer(vecLDc), 
        as.integer(cusumVecLG), as.integer(cusumVecLD), as.integer(cusumVecLK), as.integer(cusumVecLGc),
        as.integer(cusumVecLDc),        
        as.double(c(blockSizeProbG, blockSizeProbD, blockSizeProbK, blockSizeProbGc, blockSizeProbDc)), 
        as.integer(c(maxBSG,maxBSD,maxBSK,maxBSGc,maxBSDc)),                        
        as.double(tuneCalpha), as.double(tuneSigma2), as.double(tuneCbeta), as.double(tuneAlpha), 
        as.double(tuneSigma2R), 
        as.double(tuneR), as.double(tuneCpsi), as.double(tuneCbCor), as.double(tuneOmega), 
        as.double(tuneComega),                                    
        as.double(pimu), as.double(cetaParams), as.double(pisigma), as.integer(HNca), 
        as.double(calphaParams),            
        as.integer(HNsg), as.double(sigmaParams), as.double(piphi), as.integer(HNcpsi), 
        as.double(cpsiParams),
        as.double(cetaCorParams),as.integer(HNco),as.double(comegaParams),as.double(c(pinu,pifi)),
        as.integer(HNscor),as.double(sigmaCorParams),
        as.double(c(tau,tol)), as.integer(FT), as.double(deviance),as.integer(isDz2),            
        as.integer(c(cont,LASTsw,LASTWB)),as.double(LASTAll))}
    if (mcm==6){
        out<-.C("longmultg",
        as.integer(seed), as.character(StorageDir), 
        as.integer(c(sweeps,burn,thin,n,p,N,MVLD,mcm)),
        as.integer(niVec), as.integer(cusumniVec), as.integer(intime2), as.integer(intime), 
        as.integer(c(niMax,LUT)), as.integer(FUT),
        as.integer(cusumC), as.double(C), as.double(c(t(Y))), as.double(t(X)), as.double(Z), 
        as.double(Xc),as.double(Zc),
        as.integer(c(LG,LD,LK,LGc,LDc)), as.integer(c(NG,ND2,NK,NGc,NDc)), 
        as.integer(vecLG), as.integer(vecLD), as.integer(vecLK), as.integer(vecLGc),as.integer(vecLDc), 
        as.integer(cusumVecLG), as.integer(cusumVecLD), as.integer(cusumVecLK), as.integer(cusumVecLGc),
        as.integer(cusumVecLDc),
        as.double(c(blockSizeProbG, blockSizeProbD, blockSizeProbK, blockSizeProbGc, blockSizeProbDc)), 
        as.integer(c(maxBSG,maxBSD,maxBSK,maxBSGc,maxBSDc)),                        
        as.double(tuneCalpha), as.double(tuneSigma2), as.double(tuneCbeta), as.double(tuneAlpha), 
        as.double(tuneSigma2R), 
        as.double(tuneR), as.double(tuneCpsi), as.double(tuneCbCor), as.double(tuneOmega), 
        as.double(tuneComega),  
        as.double(pimu), as.double(cetaParams), as.double(pisigma), as.integer(HNca), 
        as.double(calphaParams),            
        as.integer(HNsg), as.double(sigmaParams), as.double(piphi), as.integer(HNcpsi), 
        as.double(cpsiParams),
        as.double(cetaCorParams),as.integer(HNco),as.double(comegaParams),as.double(c(pinu,pifi)),
        as.integer(HNscor),as.double(sigmaCorParams),
        as.double(tau), as.integer(FT), as.double(deviance),as.integer(isDz2),            
        as.integer(c(cont,LASTsw,LASTWB)),as.double(LASTAll),
        as.integer(H), as.double(DPparams))}   
    if (mcm==7){
        out<-.C("longmultgv",
        as.integer(seed), as.character(StorageDir), 
        as.integer(c(sweeps,burn,thin,n,p,N,MVLD,mcm)),
        as.integer(niVec), as.integer(cusumniVec), as.integer(intime2), as.integer(intime), 
        as.integer(c(niMax,LUT)), as.integer(FUT),
        as.integer(cusumC), as.double(C), as.double(c(t(Y))), as.double(t(X)), as.double(Z), 
        as.double(Xc),as.double(Zc),
        as.integer(c(LG,LD,LK,LGc,LDc)), as.integer(c(NG,ND2,NK,NGc,NDc)), 
        as.integer(vecLG), as.integer(vecLD), as.integer(vecLK), as.integer(vecLGc),as.integer(vecLDc), 
        as.integer(cusumVecLG), as.integer(cusumVecLD), as.integer(cusumVecLK), as.integer(cusumVecLGc),as.integer(cusumVecLDc),          
        as.double(c(blockSizeProbG, blockSizeProbD, blockSizeProbK, blockSizeProbGc, blockSizeProbDc)), 
        as.integer(c(maxBSG,maxBSD,maxBSK,maxBSGc,maxBSDc)),                        
        as.double(tuneCalpha), as.double(tuneSigma2), as.double(tuneCbeta), as.double(tuneAlpha), as.double(tuneSigma2R), 
        as.double(tuneR), as.double(tuneCpsi), as.double(tuneCbCor), as.double(tuneOmega), as.double(tuneComega),  
        as.double(pimu), as.double(cetaParams), as.double(pisigma), as.integer(HNca), as.double(calphaParams),            
        as.integer(HNsg), as.double(sigmaParams), as.double(piphi), as.integer(HNcpsi), as.double(cpsiParams),
        as.double(cetaCorParams),as.integer(HNco),as.double(comegaParams),as.double(c(pinu,pifi)),
        as.integer(HNscor),as.double(sigmaCorParams),
        as.double(tau), as.integer(FT), as.double(deviance),as.integer(isDz2),            
        as.integer(c(cont,LASTsw,LASTWB)), as.double(LASTAll),
        as.integer(G), as.double(DPparams))}   
    #Output  
    if (!is.null(C)) C<-t(C)
    loc1<-31 
    loc2<-60 
    fit <- list(call=call,call2=call2,formula=formula,seed=seed,p=p,d=p*(p-1)/2,
                data=data,lag=lag,Y=Y,
                X=X,Xknots=Xknots,LG=LG,NG=NG,
                Z=Z,Zknots=Zknots,LD=LD,ND=ND,ND2=ND2,
                storeMeanVectorX=storeMeanVectorX,
                storeMeanVectorZ=storeMeanVectorZ,
                storeIndicatorX=storeIndicatorX,
                storeIndicatorZ=storeIndicatorZ,
                assignX=assignX,
                assignZ=assignZ,
                labelsX=labelsX,
                labelsZ=labelsZ,
                countX=countX,
                countZ=countZ,
                varsY=varsY,                
                varsX=varsX,
                varsZ=varsZ,
                is.Dx=is.Dx,
                is.Dz=is.Dz,
                which.SpecX=which.SpecX,
                which.SpecZ=which.SpecZ,
                formula.termsX=formula.termsX,
                formula.termsZ=formula.termsZ,                
                nSamples=nSamples,
                totalSweeps=sweeps,
                mcpar=c(as.integer(burn+1),as.integer(seq(from=burn+1,by=thin,length.out=nSamples)[nSamples]),as.integer(thin)),                                
                mcm=mcm,H=H,G=G,                
                tuneCalpha=c(tuneCalpha,out[[loc1+0]][1:p]),        
                tuneSigma2=c(tuneSigma2,out[[loc1+1]][1:p]),
                tuneCbeta=c(tuneCbeta,out[[loc1+2]][1]),                                
                tuneAlpha=c(tuneAlpha,out[[loc1+3]][1:(p*ND2)]),                
                tuneSigma2R=c(tuneSigma2R,out[[loc1+4]][1]),
                tuneR=c(tuneR,out[[loc1+5]][1:LUT]),
                tuneCpsi=c(tuneCpsi,out[[loc1+6]][1:(p*p)]),                
	            tuneCbCor=c(tuneCbCor,out[[loc1+7]][1]),	            
	            tuneOmega=c(tuneOmega,out[[loc1+8]][1:NDc]),	            
                tuneComega=c(tuneComega,out[[loc1+9]][1]),
                deviance=c(out[[loc2]][1:2]),
                nullDeviance=nullDeviance,
                DIR=StorageDir,  
                out=out,
                LUT=LUT,
                SUT=SUT,
                C=C,Cknots=Cknots,LK=LK,NK=NK,
                Xc=Xc,Xcknots=Xcknots,LGc=LGc,NGc=NGc,
                Zc=Zc,Zcknots=Zcknots,LDc=LDc,NDc=NDc,                
                storeMeanVectorC=storeMeanVectorC,
                storeMeanVectorXc=storeMeanVectorXc,
                storeMeanVectorZc=storeMeanVectorZc,                               
                storeIndicatorC=storeIndicatorC,
                storeIndicatorXc=storeIndicatorXc,
                storeIndicatorZc=storeIndicatorZc, 
                assignC=assignC,          
                assignXc=assignXc,
                assignZc=assignZc,
                labelsC=labelsC,                
                labelsXc=labelsXc,
                labelsZc=labelsZc,                                               
                countC=countC,                            
                countXc=countXc,
                countZc=countZc,
                varsC=varsC,                
                varsXc=varsXc,
                varsZc=varsZc,
                varTime=varTime,  
                is.Dc=is.Dc,
                is.Dxc=is.Dxc,
                is.Dzc=is.Dzc,  
                which.SpecC=which.SpecC,                
                which.SpecXc=which.SpecXc,
                which.SpecZc=which.SpecZc, 
                formula.termsC=formula.termsC,
                formula.termsXc=formula.termsXc,
                formula.termsZc=formula.termsZc,                 
                HNca=HNca,
                HNcpsi=HNcpsi,                
                HNsg=HNsg,
                HNco=HNco,
                HNscor=HNscor,
                niVec=niVec,
                intime=intime,
                NC=0,
                FT=FT,
                qCont=0)                                                    
    class(fit) <- 'mvrm'
    return(fit)
} 

#psiAll <- mvrm2mcmc(x,"psi")
#alphaAll <- mvrm2mcmc(x,"alpha") 
#sigma2All <- mvrm2mcmc(x,"sigma2")
#RAll <- mvrm2mcmc(x,"R")    

postSigma<-function(x, time, psi, alpha, sigma2, R, samples, ...){
    p<-x$p
    ifelse(missing(time), t<-x$SUT, t<-time)
    T<-length(t)
    postS<-matrix(0,nrow=T*p,ncol=T*p)    

    # Vars
    varsd<-as.list(substitute(list(...)))[-1]
    if (x$varTime %in% names(varsd)) 
        varsd <-varsd[-which(names(varsd)==x$varTime)]
    if ("lag" %in% names(varsd)) 
        varsd <-varsd[-which(names(varsd)=="lag")]

    comp.des.mats<-1
    if (comp.des.mats==1){
		comp.des.mats<-0
        # Zmat
        vars<-x$varsZ
        wv<-NULL 
        newdata<-t
        for (v in 1:length(vars)){
	    	if (vars[[v]] %in% names(varsd)){
                min1<-eval(varsd[[which(vars[[v]]==names(varsd))]])
                newdata<-cbind(newdata,rep(min1,T))
                wv<-c(wv,v)
		    }
        }
        if (! (length(wv)+1) == length(vars)) stop("insufficient input on covariates")
        newdata <- as.data.frame(newdata) 
        colnames(newdata)<-c(x$varTime,vars[wv])
        terms.reform<-NULL
        k<-0
        for (i in 1:length(x$formula.termsZ)){
            term<-x$formula.termsZ[i]
            if (!i %in%  which(unlist(x$which.SpecZ)==-99)){
                k<-k+1
                if (!grepl("knots",term)){
                    term<-substr(term,1,nchar(term)-1)
     		        term<-paste(term,",knots= knots[[",k,"]])")
			    }
	        }
            terms.reform<-c(terms.reform,term)
        }
        formula2<-reformulate(terms.reform)
        if (length(x$data)>0){
            nd<-x$data[0,match(colnames(newdata),colnames(x$data)),drop=FALSE]
            for (j in 1:dim(nd)[2])
                nd[,j]<-drop(nd[,j])  
            nd[1:NROW(newdata),] <- newdata
        }else{nd<-newdata}
        Zmat<-DM(formula=formula2,data=nd,n=NROW(nd),knots=x$Zknots,meanVector=x$storeMeanVectorZ,indicator=x$storeIndicatorZ,centre=TRUE)$X
        Zmat<-as.matrix(Zmat)[,-1]    
        # Cmat  
        vars<-x$varsC
        wv<-NULL 
        C<-cbind(rep(seq(2,T,1),seq(1,T-1,1)), unlist(sapply(1:(T-1), function(i) seq(1,i,1))))
        lag<-t[C[,1]]-t[C[,2]]
        newdata<-lag
        for (v in 1:length(vars)){
	    	if (vars[[v]] %in% names(varsd)){
                min1<-eval(varsd[[which(vars[[v]]==names(varsd))]])
                newdata<-cbind(newdata,rep(min1,T))
                wv<-c(wv,v)
		    }
        }
        if (! (length(wv)+1) == length(vars)) stop("insufficient input on covariates")
        newdata<-as.data.frame(newdata)
        colnames(newdata)<-c("lag",vars[wv])
        newdata <- as.data.frame(newdata) 
        terms.reform<-NULL
        k<-0
        for (i in 1:length(x$formula.termsC)){
            term<-x$formula.termsC[i]
            if (!i %in%  which(unlist(x$which.SpecC)==-99)){
                k<-k+1
                if (!grepl("knots",term)){
                    term<-substr(term,1,nchar(term)-1)
     		        term<-paste(term,",knots= knots[[",k,"]])")
			    }
	        }
            terms.reform<-c(terms.reform,term)
        }
        formula2<-reformulate(terms.reform)
        if (length(x$data)>0){
	    	dataNew<-as.data.frame(cbind(x$data,lag=rnorm(dim(x$data)[1])))
            nd<-dataNew[0,match(colnames(newdata),colnames(dataNew)),drop=FALSE]
            for (j in 1:dim(nd)[2])
                nd[,j]<-drop(nd[,j])
            nd[1:NROW(newdata),] <- newdata
        }else{nd<-newdata}    
        Cmat<-DM(formula=formula2,data=nd,n=NROW(newdata),knots=x$Cknots,meanVector=x$storeMeanVectorC,indicator=x$storeIndicatorC,centre=TRUE)$X   
        Cmat<-as.matrix(Cmat)
    }
    
    # #
    Sij<-matrix(0,p,p)
    RA<-diag(rep(1,p))
    Li<-diag(rep(1,p*T))
    Di<-matrix(0,p*T,p*T)
    
    if (missing(samples)) samples <- 1:x$nSamples
    
    for (sw in samples){		        
	    buildPhi<-NULL
		for (r1 in 1:p){
	        for (r2 in 1:p){
		        Pair<-(r1-1)*p + r2
		        psiA<-psi[sw,(1+((Pair-1)*x$LK)):(Pair*x$LK)]		        
		        buildPhi<-cbind(buildPhi,Cmat%*%matrix(psiA))
            }
		}				        
		move<-1
        if (T > 1){
            for (j in 2:T){
                for (k in 1:(j-1)){
                    Li[(j*p-p+1):(j*p),(k*p-p+1):(k*p)] <- -matrix(buildPhi[move,],p,p,byrow=TRUE) 		
                    move <- move + 1
			    }
		    }
		}
		Sr<-NULL
		for (r in 1:p){
			s2<-sigma2[sw,r]
            alphaA<-alpha[sw,(1+((r-1)*x$LD)):(r*x$LD)]                                                            
            Sr<-cbind(Sr,sqrt(s2*exp(Zmat%*%matrix(alphaA))))
        }
        for (j in 1:T){
            diag(Sij) <- Sr[j,]  
            pick.time<-which(t[j]==x$SUT)#j ; x$intime[cusumniVec[subject]+j]+1
            RA[lower.tri(RA)] <- R[sw,(1+(pick.time-1)*p*(p-1)/2):(pick.time*p*(p-1)/2)]
            RA[upper.tri(RA)] <- t(RA)[upper.tri(RA)] 
            Di[((j-1)*p+1):(j*p),((j-1)*p+1):(j*p)] <- Sij %*% RA %*% Sij
        }       
        SigmaSW <- solve(Li)%*%Di%*%solve(t(Li))
        postS <- postS + SigmaSW
	}
	return(postS/length(samples))
}

postSigma2<-function(x, subject, psi, alpha, sigma2, R, samples, ...){
    niVec<-x$niVec
    p<-x$p
    T<-niVec[subject]
    postS<-matrix(0,nrow=T*p,ncol=T*p)    
    cusumniVec<-c(0,cumsum(niVec))
    cusumC<-c(0,cumsum(niVec*(niVec-1)/2))                
    Cmat<-x$C[(cusumC[subject]+1):cusumC[subject+1],]      
    Zmat<-x$Z[(cusumniVec[subject]+1):cusumniVec[subject+1],-1]        
    Sij<-matrix(0,p,p)
    RA<-diag(rep(1,p))
    Li<-diag(rep(1,p*T))
    Di<-matrix(0,p*T,p*T)
    if (missing(samples)) samples <- 1:x$nSamples
    for (sw in samples){		        
	    buildPhi<-NULL
		for (r1 in 1:p){
	        for (r2 in 1:p){
		        Pair<-(r1-1)*p + r2
		        psiA<-psi[sw,(1+((Pair-1)*x$LK)):(Pair*x$LK)]		        
		        buildPhi<-cbind(buildPhi,Cmat%*%matrix(psiA))
            }
		}				        
		move<-1
        for (j in 2:T){
            for (k in 1:(j-1)){
                Li[(j*p-p+1):(j*p),(k*p-p+1):(k*p)] <- -matrix(buildPhi[move,],p,p,byrow=TRUE) 		
                move <- move + 1
			}
		}
		Sr<-NULL
		for (r in 1:p){
			s2<-sigma2[sw,r]
            alphaA<-alpha[sw,(1+((r-1)*x$LD)):(r*x$LD)]                                                            
            Sr<-cbind(Sr,sqrt(s2*exp(Zmat%*%matrix(alphaA))))
        }
        for (j in 1:T){
            diag(Sij) <- Sr[j,]  
            pick.time<-x$intime[cusumniVec[subject]+j]+1              
            RA[lower.tri(RA)] <- R[sw,(1+(pick.time-1)*p*(p-1)/2):(pick.time*p*(p-1)/2)]
            RA[upper.tri(RA)] <- t(RA)[upper.tri(RA)]                                    
            Di[((j-1)*p+1):(j*p),((j-1)*p+1):(j*p)] <- Sij %*% RA %*% Sij
        }       
        SigmaSW <- solve(Li)%*%Di%*%solve(t(Li))
        postS <- postS + SigmaSW
	}
	return(postS/length(samples))
}

plot3.bcmg<-function(x, model="mean", centre=mean, quantiles=c(0.10,0.90), plotEmptyCluster=FALSE, plotOptions=list(), ...){
    t1<-unique(x$Xc[,2]) # this line is not needed
    t1<-x$SUT
    if (model=="mean"){
        meanReg<-array(0,dim=c(x$nSamples,x$H,x$LUT)) 
        uXc<-unique(x$Xc)
        file1 <- paste(x$DIR,"BNSP.eta.ls.txt",sep="")
        file2 <- paste(x$DIR,"BNSP.clusters.txt",sep="")
        if (!file.exists(file1) && !file.exists(file2)) quiet(ls.mvrm(x))
        eta<-array(unlist(read.table("BNSP.eta.ls.txt")),dim=c(x$nSamples,x$H,(x$LGc+1)))
        tabs<-table(array(unlist(read.table("BNSP.clusters.txt"))))
        labs<-array(0,x$H)
        labs[as.numeric(names(tabs))]<-tabs
        for (i in 1:x$nSamples)
            meanReg[i,,] <- eta[i,,]%*%t(uXc)
        if (x$FT==1) meanReg <- tanh(meanReg)
        centreM<- apply(meanReg,c(2,3),centre)
        QM<-NULL; SQM<-matrix(NA,ncol=2,nrow=NROW(centreM))
        if (!is.null(quantiles)){ 
            QM<-apply(meanReg,c(2,3),quantile,probs=quantiles,na.rm=TRUE)
            SQM<-cbind(QLM=c(t(QM[1,,])),QUM=c(t(QM[2,,])))
	    }
        dataM<-data.frame(group=factor(rep(1:x$H,each=x$LUT)),size=rep(labs,each=x$LUT),t=rep(t1,x$H),centre=c(t(centreM)),QLM=SQM[,1],QUM=SQM[,2])		
        dataM<-dataM[with(dataM, order(-size, t)),]
        #dataM$group<-factor(rep(seq(1,x$H,1),each=x$LUT))# not needed
        labs2<-unique(dataM$size)
        if (!plotEmptyCluster) dataM<-dataM[dataM$size > 0,]
        plotElM<-geom_line(aes_string(x="t",y="centre", group="group",col="group"),alpha=1.0,cex=1.3)
        if (!is.null(quantiles))  
            plotElM<-c(geom_ribbon(aes_string(x="t",ymin="QLM",ymax="QUM",group="group",fill="group",col="group"),alpha =0.2),plotElM)
        ggM<-ggplot(data=dataM)
        plotM <- ggM + plotElM + ylab(expression(mu[t])) + guides(fill=FALSE) + 
                 scale_color_discrete(name = "group [size]", labels = paste("[",labs2,"]",sep="")) + plotOptions  
    }
    if (model=="stdev"){
	    fitV<-matrix(0,nrow=x$nSamples,ncol=length(t1))
	    uZ<-unique(x$Z)
		alphaFN <- file.path(paste(x$DIR,"BNSP.alpha.txt",sep="")) 
        aFile<-file(alphaFN,open="r")
		sigma2FN <- file.path(paste(x$DIR,"BNSP.sigma2.txt",sep="")) 
        s2File<-file(sigma2FN,open="r")
        for (i in 1:x$nSamples){
        	alpha<-scan(aFile,what=numeric(),n=x$LD,quiet=T)
		    s2<-scan(s2File,what=numeric(),n=1,quiet=T)
            fitV[i,]<-sqrt(s2*exp(uZ[,-1]%*%matrix(alpha)))
		}
		close(aFile)
		close(s2File)	
		centreM<-apply(fitV,2,centre)
        QM<-NULL; SQM<-matrix(NA,ncol=2,nrow=NROW(centreM))
        if (!is.null(quantiles)){ 
            QM<-apply(fitV,2,quantile,probs=quantiles,na.rm=TRUE)
            SQM<-cbind(QLM=QM[1,],QUM=QM[2,])
	    }
        dataM<-data.frame(t=t1,centre=centreM,QLM=SQM[,1],QUM=SQM[,2])		
        plotElM<-geom_line(aes_string(x="t",y="centre"),alpha=1.0,cex=1.3)
        if (!is.null(quantiles))  
            plotElM<-c(geom_ribbon(aes_string(x="t",ymin="QLM",ymax="QUM"),alpha =0.2),plotElM)
        ggM<-ggplot(data=dataM)
        plotM <- ggM + plotElM + ylab(expression(sigma[t])) + plotOptions  
	}
	return(plotM)
}

Try the BNSP package in your browser

Any scripts or data that you put into this service are public.

BNSP documentation built on May 31, 2023, 7:05 p.m.