#
# mining.r (Coal-Mining Disasters Example)
#

# set constants, parameters and data
m <- 100
t <- 15
n <- 112
y <- c(4, 5, 4, 1, 0, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6, 3, 3,
    5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5, 2, 2, 3, 4, 2,
    1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0, 1, 0, 1, 1, 0, 0, 3, 1,
    0, 3, 2, 2, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0,
    0, 1, 1, 0, 2, 3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0,
    0, 1, 4, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1)
a1 <- 0.5
a2 <- 0.5
c1 <- 0
c2 <- 0
d1 <- 1
d2 <- 1
# calculate partial sums that are frequently needed
partsum <- (0:n)
for (i in 2:(n+1))
{
    partsum[i] <- partsum[i-1]+y[i-1]
}
# set initial data (arbitrary choices)
k <- (1:m)
l1 <- (1:m)
l2 <- (1:m)
b1 <- (1:m)
b2 <- (1:m)
for (i in 1:m)
{
    b1[i] <- 0.5
    b2[i] <- 0.5
    k[i] <- 56
}
# execute the algorithm t times
for (j in 1:t)
{
    # create m independent sequences
    for (l in 1:m)
    {
        # step for lambda_1
        l1[l] <- rgamma(1, a1+partsum[k[l]+1],
            b1[l]/(k[l]*b1[l]+1))
        # step for lambda_2
        l2[l] <-
            rgamma(1, a2+partsum[n+1]-partsum[k[l]+1],
                b2[l]/((n-k[l])*b2[l]+1))
        # step for b_1
        b1[l] <- 1/rgamma(1, a1+c1, d1/(l1[l]*d1+1))
        # step for b_2
        b2[l] <- 1/rgamma(1, a2+c2, d2/(l2[l]*d2+1))
        # step for k
        # calculate the inverses due to overflows in the exponentials
        invers <- (1:n)
        for (i in 1:n)
        {
            invers[i] <- 0
            for (p in 1:n)
                invers[i] <- invers[i] +
                    exp((l2[l]-l1[l])*(p-i)+
                    (partsum[p+1]-partsum[i+1])*log(l1[l]/l2[l]))
        }
        # generate vector with probabilities
        prob <- 1/invers
        # sample
        k[l] <- sample(n,1,prob)
    }
}
# Monte Carlo Integration, here only for k
marg <- (1:n)
for (i in 1:n)
{
    marg[i] <- 0
}
for (l in 1:m)
{
    # use the same scheme as in the algorithm
    invers <- (1:n)
    for (i in 1:n)
    {
        invers[i] <- 0
        for (p in 1:n)
            invers[i] <- invers[i] +
                exp((l2[l]-l1[l])*(p-i)+
                (partsum[p+1]-partsum[i+1])*log(l1[l]/l2[l]))
    }
    prob <- 1/invers
    marg <- marg + prob/m
}
# calculate expected value
print(sum((1:n)*marg))
# plot the marginal distribution
plot((1:n),marg)