Mixture of Poisson Models: Earthquake

rm(list=ls())

Problem Determination

## Data (including year 2007-2016)
dat <- read.table("earthquakes.txt", header = FALSE)
names(dat) <- c("year","earthquake")

plot(dat,type = "b", xlab = "Year", ylab = "Number",
     main = "Number of major earthquakes")

earthquake.table <- table(dat$earthquake)
plot(earthquake.table, xlim = c(0,45), ylim = c(0,30), xlab = "Number", ylab = "Frequency", main = "Frequency of numbers of major earthquakes")
points(0:45, dpois(0:45, lambda = mean(dat$earthquake)) * dim(dat)[1], pch = 19)

mean(dat$earthquake)
var(dat$earthquake)

Hence overdispersion.

acf(dat$earthquake)

The function acf computes (and by default plots) estimates of the autocovariance or autocorrelation function. Hence autocorreltaion.

Mixture Model for number of Earthquake

Poisson mixture: transform

## natural to working parameters
poisMix.pn2pw <- function(nDist, lambda, delta){
    # length(lambda) = nDist
    # length(delta) = nDist - 1
    if(sum(delta) >= 1){print("sum(delta) should be < 1")
    return()}
    if(length(lambda) != nDist){
        print("length(lambda) should be nDist")
        return()
    }
    if(length(delta) != (nDist - 1)){
        print("length(delta) should be nDist - 1")
        return()
    }
    eta <- log(lambda)
    tau <- log(delta / (1 - sum(delta)))
    return(list(eta = eta, tau = tau))
}

## Poisson mixture: transform
## working to natural parameters
poisMix.pw2pn <- function(nDist,eta,tau){
    if(nDist == 1){return(exp(eta))}
    if(length(eta) != nDist){
        print("length(lambda) should be nDist")
    return()}
    if(length(tau) != (nDist-1)){
        print("length(delta) should be nDist-1")
        return()
    }
    lambda <- exp(eta)
    delta <- exp(tau) / (1 + sum(exp(tau)))
    delta <- c(1 - sum(delta), delta)
    return(list(lambda = lambda, delta = delta))
}

## negative log likelihood function
negaLogL <- function(theta, nDist, y){
    if(nDist == 1){
        return(- sum(dpois(y, lambda = exp(theta), log = TRUE)))
    }
    eta <- theta[1: nDist]
    tau <- theta[(nDist + 1): (2 * nDist - 1)]
    naturalPars <- poisMix.pw2pn(nDist, eta, tau)
    n <- length(y)
    nll <- 0  # negative log likelihood
    for(i in 1:n){
        nll <- nll - log(sum(naturalPars$delta * 
                             dpois(y[i], lambda = naturalPars$lambda)))
    }
    return(nll)
}

Estimation with one distribution

m_HMM_mix_1 <- 1; ## No of states
## Initial values
lambda_mix.0_1 <- mean(dat$earthquake); 
delta_mix.0_1 <- c()
## Working parameters
workPars_mix_1 <- poisMix.pn2pw(m_HMM_mix_1, lambda_mix.0_1, delta_mix.0_1)
theta_mix.0_1 <- c(workPars_mix_1$eta, workPars_mix_1$tau)
## MLE
opt_mix_1 <- nlminb(theta_mix.0_1, negaLogL, nDist = m_HMM_mix_1, y = dat$earthquake)
sprintf("Max-Likelihood estimation with one distribution")
print(opt_mix_1)
## Natural parameters
naturalPars_mix_1 <- poisMix.pw2pn(m_HMM_mix_1, opt_mix_1$par[1: m_HMM_mix_1], 
                                   opt_mix_1$par[(m_HMM_mix_1 + 1) : (2 * m_HMM_mix_1 - 1)])
sprintf("Natural parameters by estimation with one distribution")
print(naturalPars_mix_1)

Estimation with two distributions

m_HMM_mix_2 <- 2; ## No of states
## Initial values
lambda_mix.0_2 <- mean(dat$earthquake) * c(1/2, 3/2);
delta_mix.0_2 <- c(0.5)
## Working parameters
workPars_mix_2 <- poisMix.pn2pw(m_HMM_mix_2, lambda_mix.0_2, delta_mix.0_2)
theta_mix.0_2 <- c(workPars_mix_2$eta, workPars_mix_2$tau)
## MLE
opt_mix_2 <- nlminb(theta_mix.0_2, negaLogL, nDist = m_HMM_mix_2, y = dat$earthquake)
sprintf("Max-Likelihood estimation with two distributions")
print(opt_mix_2)
## Natural parameters
naturalPars_mix_2 <- poisMix.pw2pn(m_HMM_mix_2, opt_mix_2$par[1: m_HMM_mix_2],
                                   opt_mix_2$par[(m_HMM_mix_2 + 1): (2 * m_HMM_mix_2 - 1)])
sprintf("Natural parameters by estimation with two distributions")
print(naturalPars_mix_2)

Estimation with 3 distributions

m_HMM_mix_3 <- 3;
lambda_mix.0_3 <- mean(dat$earthquake) * c(1/2, 1, 3/2);
delta_mix.0_3 <- c(1/3, 1/3)
workPars_mix_3 <- poisMix.pn2pw(m_HMM_mix_3, lambda_mix.0_3, delta_mix.0_3)
theta_mix.0_3 <- c(workPars_mix_3$eta, workPars_mix_3$tau)
opt_mix_3 <- nlminb(theta_mix.0_3, negaLogL, nDist = m_HMM_mix_3, y = dat$earthquake)
sprintf("Max-Likelihood estimation with three distributions")
print(opt_mix_3)
naturalPars_mix_3 <- poisMix.pw2pn(m_HMM_mix_3, opt_mix_3$par[1: m_HMM_mix_3],
                                   opt_mix_3$par[(m_HMM_mix_3 + 1): (2 * m_HMM_mix_3 - 1)])
sprintf("Natural parameters by estimation with three distributions")
print(naturalPars_mix_3)

Estimation with 4 distributions

m_HMM_mix_4 <- 4;
lambda_mix.0_4 <- mean(dat$earthquake) * c(1/4, 3/4, 5/4, 7/4);
delta_mix.0_4 <- c(1, 1, 1) / 4
workPars_mix_4 <- poisMix.pn2pw(m_HMM_mix_4, lambda_mix.0_4, delta_mix.0_4)
theta_mix.0_4 <- c(workPars_mix_4$eta, workPars_mix_4$tau)
opt_mix_4 <- nlminb(theta_mix.0_4, negaLogL, nDist = m_HMM_mix_4, y = dat$earthquake)
sprintf("Max-Likelihood estimation with four distributions")
print(opt_mix_4)
naturalPars_mix_4 <- poisMix.pw2pn(m_HMM_mix_4, opt_mix_4$par[1: m_HMM_mix_4],
                                   opt_mix_4$par[(m_HMM_mix_4 + 1): (2 * m_HMM_mix_4 - 1)])
sprintf("Natural parameters by estimation with four distributions")
print(naturalPars_mix_4)

Plot the result

plot(earthquake.table / sum(earthquake.table), xlim = c(0,45), axes = FALSE,
     xlab = "Number", ylab = "Density", 
     main = "Poisson Mixture Model for numbers of major earthquakes")
axis(1); axis(2)
x <- 0:45
mixDist <- function(x, naturalPars){
    ## Function to calculate mixture dist
    sum(naturalPars$delta * dpois(x, lambda = naturalPars$lambda))
}  
lines(x, dpois(x, lambda = naturalPars_mix_1), 
      type = "b", col = 1, pch = 19)
lines(x, sapply(x, mixDist, naturalPars = naturalPars_mix_2), 
      type = "b", col = 2, pch = 19)
lines(x, sapply(x, mixDist, naturalPars = naturalPars_mix_3), 
      type = "b", col = 3, pch = 19)
lines(x, sapply(x, mixDist, naturalPars =naturalPars_mix_4), 
      type = "b", col = 4, pch = 19)
legend("topright", col = seq(1, 4), lwd = 2, 
       c("1 distritution", "2 distritutions", 
         "3 distritutions", "4 distritutions"))

Choose the most suitable model using AIC

AIC_poisMix <- 2 * c(opt_mix_1$objective, opt_mix_2$objective,
                     opt_mix_3$objective, opt_mix_4$objective) +
               2 * c(length(opt_mix_1$par), length(opt_mix_2$par),
                     length(opt_mix_3$par), length(opt_mix_4$par))
print(AIC_poisMix)

Hence we choose 3-state model.

Choose the most suitable model using likelihood ratio test

1 - pchisq(2 * (opt_mix_1$objective - opt_mix_2$objective),
           df = length(opt_mix_2$par)-length(opt_mix_1$par))
1 - pchisq(2 * (opt_mix_2$objective - opt_mix_3$objective),
           df=length(opt_mix_3$par)-length(opt_mix_2$par))
1 - pchisq(2 * (opt_mix_3$objective - opt_mix_4$objective),
           df=length(opt_mix_4$par) - length(opt_mix_3$par))

Hence same conclusion as for AIC.

Hidden Markov Model for number of Earthquake

source("A1.R")
y <- dat$earthquake

Hidden Markov Model with 1 State

m_HMM_1 <- 1
lambda_HMM.0_1 <- mean(y)
gamma_HMM.0_1 <- 0
opt_HMM_1 <- pois.HMM.mle(y, m_HMM_1, lambda_HMM.0_1, gamma_HMM.0_1)

Hidden Markov Model with 2 State

m_HMM_2 <- 2
lambda_HMM.0_2 <- quantile(y, c(0.25, 0.75))
gamma_HMM.0_2 <- matrix(0.05, ncol = m_HMM_2, nrow = m_HMM_2)
diag(gamma_HMM.0_2) <- 1 - (m_HMM_2 - 1) * gamma_HMM.0_2[1,1]
opt_HMM_2 <- pois.HMM.mle(y, m_HMM_2, lambda_HMM.0_2, gamma_HMM.0_2)

Hidden Markov Model with 3 State

m_HMM_3 <- 3
lambda_HMM.0_3 <- quantile(y, c(0.25, 0.5, 0.75))
gamma_HMM.0_3 <- matrix(0.05, ncol = m_HMM_3, nrow = m_HMM_3)
diag(gamma_HMM.0_3) <- 1 - (m_HMM_3 - 1) * gamma_HMM.0_3[1, 1]
opt_HMM_3 <- pois.HMM.mle(y, m_HMM_3, lambda_HMM.0_3, gamma_HMM.0_3)

Hidden Markov Model with 4 State

m_HMM_4 <- 4
lambda_HMM.0_4 <- quantile(y,c(0.2, 0.4, 0.6, 0.8))
gamma_HMM.0_4 <- matrix(0.025, ncol = m_HMM_4, nrow = m_HMM_4)
diag(gamma_HMM.0_4) <- 1 - (m_HMM_4 - 1) * gamma_HMM.0_4[1, 1]
opt_HMM_4 <- pois.HMM.mle(y, m_HMM_4, lambda_HMM.0_4, gamma_HMM.0_4)

AIC

AIC_HMM <- c(opt_HMM_1$AIC, opt_HMM_2$AIC, opt_HMM_3$AIC, opt_HMM_4$AIC)
logL_HMM <- -c(opt_HMM_1$mllk, opt_HMM_2$mllk, opt_HMM_3$mllk, opt_HMM_4$mllk)
m_HMM_HMM <- c(1, 2, 3, 4)
df_HMM <- m_HMM_HMM + (m_HMM_HMM^2 - m_HMM_HMM)  # lambda + gamma
cbind(df_HMM, logL_HMM, AIC_HMM)

So we choose HMM with 3 states.

## Finding the standard errors
m_HMM <- 3
## working parameters
partsWork_HMM  <- pois.HMM.pn2pw(m_HMM, opt_HMM_3$lambda, opt_HMM_3$gamma)
mod_HMM <- nlm(pois.HMM.mllk, partsWork_HMM, x = y, m = m_HMM, hessian = TRUE)  
## Organize the result
partsWork.nlm_HMM <- mod_HMM$estimate
names(partsWork.nlm_HMM) <- c("lambda1", "lambda2", "lambda3", "tau21", "tau31", "tau12", "tau32", "tau13", "tau23")
se_HMM <- sqrt(diag(solve(mod_HMM$hessian)))

## Working pars + standard error
round(cbind(partsWork.nlm_HMM, se_HMM), digits = 2)
opt_HMM_3$gamma

Parametric bootstrap

source("A2.R")
## Initailize
n <- length(y)
k <- 100
m <- 3
lambda0 <- quantile(y, c(0.25, 0.5, 0.75))
gamma0 <- matrix(0.025, ncol = m, nrow = m)
diag(gamma0) <- 1 - (m - 1) * gamma0[1,1]
## Matrices to stores bootstap result
GAMMA <- matrix(ncol = m * m, nrow = k)
Lambda <- matrix(ncol = m, nrow = k)
Delta <- matrix(ncol = m, nrow = k)
Code <- numeric(k)
## Boot strap with two different optimizers
for(i in 1:k){
    set.seed(i)
    ## generate sample
    y.sim <- pois.HMM.generate_sample(n, m, opt_HMM_3$lambda, opt_HMM_3$gamma)
    ## fit model to sample
    opt_HMM_3.tmp <- pois.HMM.mle(y.sim, m, lambda0, gamma0)
    ## Store result
    GAMMA[i, ] <- c(opt_HMM_3.tmp$gamma[1, ],
                    opt_HMM_3.tmp$gamma[2, ],
                    opt_HMM_3.tmp$gamma[3, ])
    Lambda[i, ] <- opt_HMM_3.tmp$lambda
    Delta[i, ] <- opt_HMM_3.tmp$delta
    Code[i] <- opt_HMM_3.tmp$code
    print(c(i,Code[i]))
}

NA/Inf replaced by maximum positive value

[1] 1 1
[1] 2 1
[1] 3 1

NA/Inf replaced by maximum positive value

[1] 4 1

NA/Inf replaced by maximum positive value

[1] 5 1
[1] 6 1
[1] 7 1
[1] 8 1
[1] 9 1

NA/Inf replaced by maximum positive value

[1] 10  1
[1] 11  1
[1] 12  1
[1] 13  1
[1] 14  1
[1] 15  1
[1] 16  1
[1] 17  1
[1] 18  1
[1] 19  2
[1] 20  1
[1] 21  1

NA/Inf replaced by maximum positive value

[1] 22  1
[1] 23  1
[1] 24  1
[1] 25  1
[1] 26  1
[1] 27  1
[1] 28  1
[1] 29  1
[1] 30  5
[1] 31  1
[1] 32  1

NA/Inf replaced by maximum positive value

[1] 33  1
[1] 34  1

NA/Inf replaced by maximum positive value

[1] 35  1
[1] 36  1
[1] 37  1

NA/Inf replaced by maximum positive value

[1] 38  1
[1] 39  1
[1] 40  1

NA/Inf replaced by maximum positive value

[1] 41  1
[1] 42  1
[1] 43  1
[1] 44  1
[1] 45  1
[1] 46  1
[1] 47  1
[1] 48  1

NA/Inf replaced by maximum positive value

[1] 49  1
[1] 50  1
[1] 51  1
[1] 52  1
[1] 53  1
[1] 54  1
[1] 55  1
[1] 56  1
[1] 57  1
[1] 58  1
[1] 59  1
[1] 60  1
[1] 61  1
[1] 62  1
[1] 63  1
[1] 64  1
[1] 65  1
[1] 66  1
[1] 67  1
[1] 68  1
[1] 69  1
[1] 70  1
[1] 71  1
[1] 72  1
[1] 73  1
[1] 74  1
[1] 75  1
[1] 76  1
[1] 77  1
[1] 78  1
[1] 79  1
[1] 80  1
[1] 81  1
[1] 82  1
[1] 83  1

NA/Inf replaced by maximum positive value

[1] 84  1
[1] 85  5
[1] 86  1
[1] 87  1
[1] 88  1
[1] 89  1

NA/Inf replaced by maximum positive value

[1] 90  1
[1] 91  1
[1] 92  1
[1] 93  1

NA/Inf replaced by maximum positive value

[1] 94  1
[1] 95  1
[1] 96  1
[1] 97  5
[1] 98  1
[1] 99  1
[1] 100   1

sum(Code!=1)

[1] 4

## Same as above with different optimizer
for(i in 1:k){
    set.seed(i)
    y.sim <- pois.HMM.generate_sample(n, m,
                                      opt_HMM_3$lambda, opt_HMM_3$gamma)
    lambda0 <- quantile(y.sim,c(0.25,0.5,0.75))
    opt_HMM_3.tmp <- pois.HMM.mle.nlminb(y.sim, m,lambda0, gamma0)
    GAMMA[i, ] <- c(opt_HMM_3.tmp$gamma[1, ],
                    opt_HMM_3.tmp$gamma[2, ],
                    opt_HMM_3.tmp$gamma[3, ])
    Lambda[i, ] <- opt_HMM_3.tmp$lambda
    Delta[i, ] <- opt_HMM_3.tmp$delta
    Code[i] <- opt_HMM_3.tmp$code
    print(c(i,Code[i]))
}

[1] 1 0
[1] 2 0
[1] 3 0
[1] 4 0
[1] 5 0
[1] 6 0
[1] 7 0
[1] 8 0
[1] 9 0
[1] 10  0
[1] 11  0
[1] 12  0
[1] 13  0
[1] 14  0
[1] 15  0
[1] 16  0
[1] 17  0
[1] 18  0
[1] 19  0
[1] 20  0
[1] 21  0
[1] 22  0
[1] 23  0
[1] 24  0
[1] 25  0
[1] 26  0
[1] 27  0
[1] 28  0
[1] 29  0
[1] 30  0
[1] 31  0
[1] 32  0
[1] 33  0
[1] 34  0
[1] 35  0
[1] 36  0
[1] 37  0
[1] 38  0
[1] 39  0
[1] 40  0
[1] 41  0
[1] 42  0
[1] 43  0
[1] 44  0
[1] 45  0
[1] 46  0
[1] 47  0
[1] 48  0
[1] 49  0
[1] 50  0
[1] 51  0
[1] 52  0
[1] 53  0
[1] 54  0
[1] 55  0
[1] 56  0
[1] 57  0
[1] 58  0
[1] 59  0
[1] 60  0
[1] 61  0
[1] 62  0
[1] 63  0
[1] 64  0
[1] 65  0
[1] 66  0
[1] 67  0
[1] 68  0
[1] 69  0
[1] 70  0
[1] 71  0
[1] 72  0
[1] 73  0
[1] 74  0
[1] 75  0
[1] 76  0
[1] 77  0
[1] 78  0
[1] 79  0
[1] 80  0
[1] 81  0
[1] 82  0
[1] 83  0
[1] 84  0
[1] 85  0
[1] 86  0
[1] 87  0
[1] 88  0
[1] 89  0
[1] 90  0
[1] 91  0
[1] 92  0
[1] 93  0
[1] 94  0
[1] 95  0
[1] 96  0
[1] 97  0
[1] 98  0
[1] 99  0
[1] 100   0

sum(Code!=0)

[1] 0

## Plot the results (i.e. statistical distribution of
## estimates
par(mfrow = c(1,3))
hist(Lambda[ ,1])
rug(opt_HMM_3$lambda, lwd = 2, col = 2)

some values will be clipped

hist(Lambda[ ,2])
rug(opt_HMM_3$lambda, lwd = 2, col = 2)

some values will be clipped

hist(Lambda[ ,3])
rug(opt_HMM_3$lambda, lwd = 2, col = 2)

some values will be clipped