File: gmGS.R

package info (click to toggle)
r-cran-ctmcd 1.4.1-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 1,384 kB
  • sloc: cpp: 183; ansic: 19; makefile: 2
file content (80 lines) | stat: -rw-r--r-- 3,037 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
gmGS <-
function (tmabs, te, prior, burnin, conv_pvalue = 0, conv_freq = 10, 
    niter = 10000, sampl_method = "Unif", expmethod = "PadeRBS", 
    verbose = FALSE, combmat = NULL, sampl_func = NULL) 
{
    n = ncol(tmabs)
    gmest = matrix(0, n, n)
    for (i in 1:n) {
        for (j in setdiff(1:n, i)) {
            gmest[i, j] = rgamma(1, prior[[1]][i, j], rate = prior[[2]][i])
        }
    }
    diag(gmest) = -rowSums(gmest)
    gmestList = list()
    count = 0
    if(verbose==TRUE){
      pb=txtProgressBar(min=0,max=burnin+niter,title="Burn-in",style=3)
    }
    repeat {
        if (sampl_method == "Unif") {
            draw = rNijTRiT_Unif(tmabs, te, gmest, expm(gmest, 
                method = expmethod))
        }
        else if (sampl_method == "ModRej") {
            draw = rNijTRiT_ModRej(tmabs, te, gmest)
        }
        else if (sampl_method == "Comb") {
            Mmat = (combmat == "U") * tmabs
            Umat = (combmat == "M") * tmabs
            Mdraw = rNijTRiT_ModRej(Mmat, te, gmest)
            Udraw = rNijTRiT_Unif(Umat, te, gmest, expm(gmest, 
                method = expmethod))
            draw = list()
            draw$NijT = Mdraw$NijT + Udraw$NijT
            draw$RiT = Mdraw$RiT + Udraw$RiT
        }
        else if (sampl_method == "Ext") {
            draw = sampl_func(tmabs, te, gmest)
        }
        for (i in 1:n) {
            for (j in setdiff(1:n, i)) {
                gmest[i, j] = rgamma(1, draw$NijT[i, j] + prior[[1]][i, 
                  j], rate = (draw$RiT[i] + prior[[2]][i]))
            }
        }
        diag(gmest) = 0
        diag(gmest) = -rowSums(gmest)
        count = count + 1
        if (count > burnin) {
            gmestList[[length(gmestList) + 1]] = gmest
        }
        if (verbose == TRUE) setTxtProgressBar(pb,count)
        if (conv_pvalue>0 && count > burnin + 10 && (count - burnin)%%round(niter/conv_freq) == 
            0) {
            parmat = t(matrix(unlist(gmestList), n^2))
            parmat2 = parmat[, which(colSums(parmat) > 0)]
            mcmcobj = as.mcmc(parmat2)
            heidel = heidel.diag(mcmcobj, pvalue = conv_pvalue)
            if (verbose == TRUE) {
                print("Check Convergence")
            }
            if (sum(heidel[, 1]) == length(heidel[, 1]) & sum(heidel[, 
                4]) == length(heidel[, 4])) {
                if (verbose == TRUE) {
                  print("Converged according to Heidelberger and Welch Criterion")
                }
                break
            }
            if (verbose == TRUE) {
                print("Not Converged according to Heidelberger and Welch Criterion")
            }
        }
        if (count == (burnin + niter)) break
    }
    gmall = Reduce("+", gmestList)/(count - burnin)
    colnames(gmall) = colnames(tmabs)
    rownames(gmall) = rownames(tmabs)
    est = list(par = gmall, burnin = burnin, niter = count, draws = gmestList)
    return(est)
}