Nothing
# This file contains the procedure to fit the mixed Erlang - Pareto splicing model
# to interval censored data as described in Reynkens et al. (2017) and
# Section 4.3 in Albrecher et al. (2017).
###########################################################################
# Numerical tolerance
numtol <- .Machine$double.eps^0.5
# Fit splicing model of mixed Erlang and Pareto to interval censored data
.SpliceFiticPareto <- function(L, U, censored, tsplice, M = 3, s = 1:10, trunclower = 0, truncupper = Inf, ncores = NULL,
criterium = c("BIC", "AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf, cpp = FALSE) {
tsplice <- as.numeric(tsplice)
# Repeat censored if a single number
if (length(censored) == 1) {
censored <- rep(censored, length(L))
}
######
# Check input
.spliceEM_checkInput(L=L, U=U, censored=censored, trunclower=trunclower, tsplice=tsplice,
truncupper=truncupper, eps=eps, beta_tol=beta_tol, maxiter=maxiter)
# Check input for ncores
if (is.null(ncores)) ncores <- max(detectCores()-1, 1)
if (is.na(ncores)) ncores <- 1
# Match input argument for criterium
criterium <- match.arg(criterium)
##
# Initial value for pi
# Faster to already compute it here
if (requireNamespace("interval", quietly = TRUE)) {
pi_initial <- 1-.Turnbull2(tsplice, L=L, R=U, censored=censored, trunclower=trunclower, truncupper=truncupper)$surv
} else {
pi_initial <- 1-Turnbull(tsplice, L=L, R=U, censored=censored, trunclower=trunclower, truncupper=truncupper)$surv
}
## Fit model
fit <- .spliceEM_tune(lower=L, upper=U, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
M=M, s=s, pi_initial=pi_initial, nCores=ncores, criterium=criterium, reduceM=reduceM,
eps=eps, beta_tol=beta_tol, maxiter=maxiter, cpp=cpp)
# Output as MEfit object
MEfit <- .MEoutput(fit)
# EVT part
EVTfit <- EVTfit(gamma=fit$best_model$gamma, endpoint=truncupper, sigma=NULL)
# Return correct type
type <- c("ME", ifelse(is.finite(truncupper), "TPa", "Pa"))
# Log-likelihood
loglik <- fit$best_model$loglikelihood
# Compute ICs, parameters: alpha (1 less because sums to 1), r, theta, gamma's, pi's (fit$best_model$pi has length 1 here!)
df <- 2*MEfit$M-1 + 1 + length(EVTfit$gamma) + length(fit$best_model$pi)
aic <- -2*loglik + 2*df
bic <- -2*loglik + log(length(L))*df
ic <- c(aic, bic)
names(ic) <- c("AIC", "BIC")
# Return SpliceFit object
return( SpliceFit(const=fit$best_model$pi, trunclower=trunclower, t=tsplice, type=type,
MEfit=MEfit, EVTfit=EVTfit, loglik=loglik, IC=ic) )
}
# Check input arguments for SpliceFiticPareto
.spliceEM_checkInput <- function(L, U, censored = NULL, trunclower, tsplice, truncupper, eps, beta_tol, maxiter) {
##
# Check lengths
if (length(L) != length(U)) {
stop("L and U should have equal length.")
}
if (length(trunclower) != 1) {
stop("trunclower should have length 1.")
}
if (length(tsplice) != 1) {
stop("tsplice should have length 1.")
}
if (length(truncupper) != 1 ) {
stop("truncupper should have length 1.")
}
##
# Check censored
if (!is.null(censored)) {
if (!is.logical(censored)) {
stop("censored should be a logical vector.")
}
if (length(censored) != length(L)) {
stop("censored should have the same length as L.")
}
if (sum(censored) == length(L)) {
stop("Not all observations can be censored.")
}
if (any(abs(U-L)[!censored] > numtol)) {
stop("Uncensored observations should have L=U.")
}
if (any(abs(U-L)[censored] < numtol)) {
stop("Censored observations cannot have L=U.")
}
}
##
# Check data types
if (any(!is.numeric(L))) {
stop("L should consist of numerics.")
}
if (any(!is.numeric(U))) {
stop("U should consist of numerics.")
}
if (!is.numeric(trunclower)) {
stop("trunclower should be numeric.")
}
if (!is.numeric(tsplice)) {
stop("tsplice should be numeric.")
}
if (!is.numeric(truncupper)) {
stop("truncupper should be numeric.")
}
##
# Check inequalities
if (trunclower < 0) {
stop("trunclower cannot be strictly negative.")
}
if (any(L > U)) {
stop("all elements of L should be smaller than (or equal to) the corresponding elements of U.")
}
if (trunclower > tsplice) {
stop("trunclower should be smaller than (or equal to) tsplice.")
}
if (truncupper <= tsplice) {
stop("truncupper should be strictly larger than tsplice.")
}
if (!is.finite(trunclower) ) {
stop("trunclower has to be finite.")
}
if (!is.finite(tsplice)) {
stop("tsplice has to be finite.")
}
if (any(trunclower > L)) {
stop("trunclower should be smaller than (or equal to) all elements of L.")
}
if (any(U > truncupper)) {
stop("truncupper should be larger than (or equal to) all elements of U.")
}
##
# Check if enough data in all categories
# t^l <= x_i=l_i=u_i <= t <T
ind1 <- which(!censored & (U <= tsplice))
# t^l < t < x_i=l_i=u_i <=T
ind2 <- which(!censored & (U > tsplice))
# t^l <= l_i < u_i <= t <T
ind3 <- which(censored & U <= tsplice)
# t^l < t <= l_i < u_i <=T
ind4 <- which(censored & L >= tsplice)
# t^l <= l_i < t < u_i <=T
#ind5 <- which(censored & L<tsplice & U > tsplice)
if (length(ind1) == 0 & length(ind3) == 0) {
stop("No data is given for the ME part.")
}
if (length(ind2) == 0 & length(ind4) == 0) {
stop("No data is given for the Pareto part.")
}
##
# Check input for eps
if (!is.numeric(eps) | length(eps) > 1) stop("eps should be a numeric of length 1.")
if (eps <= 0) stop("eps should be strictly positive.")
# Check input for beta_tol
if (!is.numeric(beta_tol) | length(beta_tol) > 1) stop("beta_tol should be a numeric of length 1.")
if (beta_tol > 1 | beta_tol <= 0) stop("beta_tol should be in (0,1].")
# Check input for maxiter
if (!is.numeric(maxiter) | length(maxiter) > 1) stop("maxiter should be a numeric of length 1.")
if (maxiter < 1) stop("maxiter should be at least 1.")
}
## Tune the initialising parameters M and s using a grid search over the supplied parameter ranges
.spliceEM_tune <- function(lower, upper = lower, censored = NULL, trunclower = 0, tsplice, truncupper = Inf, M = 10, s = 1, pi_initial = NULL,
nCores = detectCores(), criterium = "BIC", reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf, cpp = FALSE) {
# Censoring indicator for each observation
if (is.null(censored)) {
censored <- abs(lower-upper) > numtol & !is.na(lower != upper)
}
# All combinations of M and s
tuning_parameters <- expand.grid(M, s)
if (nCores == 1) {
# Special case when nCores=1 without any cluster
i <- 1
all_model <- foreach(i = 1:nrow(tuning_parameters), .errorhandling = 'pass') %do% {
suppressWarnings(.spliceEM_fit(lower=lower, upper=upper, censored=censored,
trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
M = tuning_parameters[i, 1], s = tuning_parameters[i, 2], pi_initial=pi_initial,
criterium=criterium, reduceM=reduceM, eps=eps, beta_tol=beta_tol, maxiter=maxiter, cpp=cpp))
}
} else {
#######
# Original code
cl <- makePSOCKcluster(nCores)
registerDoParallel(cl)
i <- 1
all_model <- foreach(i = 1:nrow(tuning_parameters), .errorhandling = 'pass') %dopar% {
suppressWarnings(.spliceEM_fit(lower=lower, upper=upper, censored=censored,
trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
M = tuning_parameters[i, 1], s = tuning_parameters[i, 2], pi_initial=pi_initial,
criterium=criterium, reduceM=reduceM, eps=eps, beta_tol=beta_tol, maxiter=maxiter, cpp=cpp))
}
stopCluster(cl)
}
# Initial values for s and M and obtained values for information criterium and M
# Return NA when an error occured in fitting
f1 <- function(x) {
if (exists(criterium, x)) {
with(x, get(criterium))
} else {
NA
}
}
crit <- sapply(all_model, f1)
f2 <- function(x) {
if (exists("M", x)) {
x$M
} else {
NA
}
}
M_res <- sapply(all_model, f2)
performances <- data.frame(tuning_parameters[,1], tuning_parameters[,2], crit, M_res)
colnames(performances) = c('M_initial', 's', criterium, 'M')
# Select model with lowest IC
best_index <- which(crit == min(crit, na.rm=TRUE))[1]
best_model <- all_model[[best_index]]
return(list(best_model = best_model, performances = performances, all_model = all_model))
}
# Fit splicing model using EM algorithm including shape adjustment and reduction of number of Erlangs
.spliceEM_fit <- function(lower, upper = lower, censored, trunclower = 0, tsplice, truncupper = Inf, M = 10, s = 1, pi_initial = NULL, criterium="AIC", reduceM = TRUE,
eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf, cpp = FALSE) {
# Get initial values
initial <- .spliceEM_initial(lower=lower, upper=upper, censored=censored, trunclower=trunclower,
tsplice=tsplice, truncupper=truncupper, M=M, s=s, pi_initial=pi_initial)
if (cpp) {
# Use C++ version
# Reduction of the shape combinations
fit <- .spliceEM_shape_red_cpp(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=initial$pi, theta=initial$theta, shape=initial$shape, beta=initial$beta, gamma=initial$gamma,
criterium=criterium, improve=reduceM, eps=eps, beta_tol=beta_tol, maxiter=maxiter, adj=FALSE)
# Subsequent adjustment and reduction of the shape combinations
fit <- .spliceEM_shape_red_cpp(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=fit$pi, theta=fit$theta, shape=fit$shape, beta=fit$beta, gamma=fit$gamma,
criterium=criterium, improve=reduceM, eps=eps, beta_tol=beta_tol, maxiter=maxiter, adj=TRUE)
} else {
# Use R version
# Reduction of the shape combinations
fit <- .spliceEM_shape_red(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=initial$pi, theta=initial$theta, shape=initial$shape, beta=initial$beta, gamma=initial$gamma,
criterium=criterium, improve=reduceM, eps=eps, beta_tol=beta_tol, maxiter=maxiter, adj=FALSE)
# Subsequent adjustment and reduction of the shape combinations
fit <- .spliceEM_shape_red(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=fit$pi, theta=fit$theta, shape=fit$shape, beta=fit$beta, gamma=fit$gamma,
criterium=criterium, improve=reduceM, eps=eps, beta_tol=beta_tol, maxiter=maxiter, adj=TRUE)
}
# Return fitted parameters of splicing model
return(list(pi = fit$pi, alpha = fit$alpha, beta = fit$beta, shape = fit$shape, theta = fit$theta, gamma = fit$gamma,
loglikelihood = fit$loglikelihood, AIC=fit$AIC, BIC=fit$BIC, M = fit$M, M_initial = M, s = s))
}
# Initial values for parameters.
# Initial value for pi can be provided. When NULL (default), it will be computed using the Turnbull estimator.
.spliceEM_initial <- function(lower, upper, censored, trunclower = 0, tsplice, truncupper = Inf, M = 10, s = 1, pi_initial = NULL) {
##
# Initial value for pi
if (is.null(pi_initial)) {
# Use Turnbull estimator
if (requireNamespace("interval", quietly = TRUE)) {
pi_initial <- 1-.Turnbull2(tsplice, L=lower, R=upper, censored=censored, trunclower=trunclower, truncupper=truncupper)$surv
} else {
pi_initial <- 1-Turnbull(tsplice, L=lower, R=upper, censored=censored, trunclower=trunclower, truncupper=truncupper)$surv
}
} else {
# Use provided value
pi <- pi_initial
}
##
# Initial values for ME parameters, see Section 4.1 in Verbelen et al. (2016)
# Data for initial step: take mean for interval censored data
uc <- lower[!censored]
ic <- (lower[censored] + upper[censored]) / 2
# Right censored
ind <- which(upper[censored] == truncupper)
ic[ind] <- (lower[censored])[ind]
# Left censored
ind <- which(lower[censored] == trunclower)
ic[ind] <- (upper[censored])[ind]
# Combine
initial_data <- c(uc, ic)
# Remove initial_data equal to 0 since they cause first shape to be 0
initial_data[initial_data == 0] <- NA
initial_dataME <- initial_data[initial_data <= tsplice]
# Initial value of theta using spread factor s
theta <- max(initial_dataME, na.rm=TRUE) / s
# Initial shapes as quantiles
shape <- unique(ceiling(quantile(initial_dataME, probs = seq(0, 1, length.out = M), na.rm = TRUE) / theta))
# Initial value for alpha's
alpha <- rep(0, length(shape))
alpha[1] <- sum(initial_dataME <= shape[1]*theta)
if (length(shape) > 1) {
for (i in 2:length(shape)) {
alpha[i] <- sum(initial_dataME <= shape[i]*theta & initial_dataME > shape[i-1]*theta)
}
}
# Keep strictly positive alpha's and corresponding shapes
shape <- shape[alpha > 0]
alpha <- alpha[alpha > 0]/sum(alpha)
# alpha to beta
t_probs <- pgamma(tsplice, shape, scale=theta) - pgamma(trunclower, shape, scale=theta)
beta <- alpha * t_probs / sum(alpha*t_probs)
##
# Initial value for gamma
# Hill estimator applied to initial data
gamma <- .Hillinternal(initial_data, threshold=tsplice)
return(list(pi=pi, gamma=gamma, theta=theta, shape=shape, alpha=alpha, beta=beta))
}
## Reduction of M based on an information criterium: AIC and BIC implemented (4.3 in Verbelen et al. 2015)
.spliceEM_shape_red <- function(lower, upper, censored, trunclower = 0, tsplice, truncupper = Inf, pi, theta, shape, beta, gamma,
criterium = "AIC", improve=TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf, adj = TRUE) {
# Fit model with initial M
n <- length(lower)
if (adj) {
# Adjust shapes
fit <- .spliceEM_shape_adj(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=pi, theta=theta, shape=shape, beta=beta, gamma=gamma,
eps=eps, beta_tol=beta_tol, maxiter=maxiter)
} else {
# Do not adjust shapes => apply EM directly
fit <- .spliceEM_fit_raw(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=pi, theta=theta, shape=shape, beta=beta, gamma=gamma,
eps=eps, beta_tol=beta_tol, maxiter=maxiter)
}
loglikelihood <- fit$loglikelihood
IC <- fit[[criterium]]
shape <- fit$shape
theta <- fit$theta
beta <- fit$beta
alpha <- fit$alpha
pi <- fit$pi
gamma <- fit$gamma
M <- length(shape)
# Try to improve IC by reducing M
while (improve && length(shape) > 1) {
# Remove shape with smallest beta
min_ind <- which.min(beta)
new_shape <- shape[-min_ind]
new_beta <- beta[-min_ind]
# Ensure that new beta's sum to 1
new_beta <- new_beta/sum(new_beta)
# Fit model with smaller M using new_beta and new_shape as starting values
if (adj) {
# Adjust shapes
fit <- .spliceEM_shape_adj(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=pi, theta=theta, shape=new_shape, beta=new_beta, gamma=gamma,
eps=eps, beta_tol=beta_tol, maxiter=maxiter)
} else {
# Do not adjust shapes => apply EM directly
fit <- .spliceEM_fit_raw(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=pi, theta=theta, shape=new_shape, beta=new_beta, gamma=gamma,
eps=eps, beta_tol=beta_tol, maxiter=maxiter)
}
new_IC <- fit[[criterium]]
# Continue improving if IC is reduced
if (new_IC < IC) {
IC <- new_IC
loglikelihood <- fit$loglikelihood
shape <- fit$shape
theta <- fit$theta
beta <- fit$beta
alpha <- fit$alpha
pi <- fit$pi
gamma <- fit$gamma
M <- length(shape)
} else {
improve <- FALSE
}
}
# Compute ICs for final model, 2*length(alpha)-1+1=2*M-1 since we estimate alpha and shape (both have length M),
# but alpha sums to 1 hence 1 degree of freedom less, and we also estimate theta.
# Finally, pi and gamma are also estimated.
df <- 2*length(alpha)-1+1+1+1
aic <- -2*loglikelihood + 2*df
bic <- -2*loglikelihood + log(n)*df
return(list(M = M, pi = pi, alpha = alpha, beta = beta, shape = shape, theta = theta, gamma = gamma,
loglikelihood = loglikelihood, AIC=aic, BIC=bic))
}
# Shape adjustments (4.2 in Verbelen et al. 2015)
.spliceEM_shape_adj <- function(lower, upper, censored, trunclower = 0, tsplice, truncupper = Inf,
pi, theta, shape, beta, gamma, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf) {
# Fit model with shape equal to the initial value (and M=length(shape))
fit <- .spliceEM_fit_raw(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=pi, theta=theta, shape=shape, beta=beta, gamma=gamma,
eps=eps, beta_tol=beta_tol, maxiter=maxiter)
loglikelihood <- fit$loglikelihood
theta <- fit$theta
shape <- fit$shape
beta <- fit$beta
alpha <- fit$alpha
pi <- fit$pi
gamma <- fit$gamma
M <- length(shape)
# Before and after are the log-likelihoods used in the outer while loop
before_loglikelihood <- -Inf
after_loglikelihood <- loglikelihood
iteration <- 1
# Improve as long as log-likelihood increases significantly
while (after_loglikelihood > before_loglikelihood + eps) {
before_loglikelihood <- after_loglikelihood
# Try changing the shapes (step 1), loop over all shapes
for (i in M:1) {
improve <- TRUE
# Increase i-th shape as long as improvement and not larger than next shape
while ( improve && (i == M || ifelse(i <= length(shape), shape[i] < shape[i+1]-1, FALSE)) ) {
# Increase i-th shape by 1
new_shape <- shape
new_shape[i] <- new_shape[i]+1
fit <- .spliceEM_fit_raw(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=pi, theta=theta, shape=new_shape, beta=beta, gamma=gamma,
eps=eps, beta_tol=beta_tol, maxiter=maxiter)
new_loglikelihood <- fit$loglikelihood
if (new_loglikelihood > loglikelihood + eps) {
loglikelihood <- new_loglikelihood
pi <- fit$pi
shape <- fit$shape
theta <- fit$theta
beta <- fit$beta
alpha <- fit$alpha
gamma <- fit$gamma
# number of shape combinations might have changed after EM algorithm if beta_tol > 0
M <- length(shape)
i <- min(i, M)
} else {
improve <- FALSE
}
}
}
# Try decreasing the shapes (step 2), loop over all shapes
for (i in 1:M) {
improve <- TRUE
# Decrease i-th shape as long as improvement and larger than 1 and not smaller than previous shape
while ( improve && ( (i == 1) || ifelse(i <= length(shape), shape[i] > shape[i-1]+1, FALSE) ) && ifelse(i <= length(shape), shape[i] > 1, FALSE)) {
# Decrease i-th shape by 1
new_shape <- shape
new_shape[i] <- new_shape[i]-1
fit <- .spliceEM_fit_raw(lower=lower, upper=upper, censored=censored, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=pi, theta=theta, shape=new_shape, beta=beta, gamma=gamma,
eps=eps, beta_tol=beta_tol, maxiter=maxiter)
new_loglikelihood <- fit$loglikelihood
if (new_loglikelihood > loglikelihood + eps) {
loglikelihood <- new_loglikelihood
pi <- fit$pi
shape <- fit$shape
theta <- fit$theta
beta <- fit$beta
alpha <- fit$alpha
gamma <- fit$gamma
# number of shape combinations might have changed after EM algorithm if beta_tol > 0
M <- length(shape)
i <- min(i, M)
} else {
improve <- FALSE
}
}
}
after_loglikelihood <- loglikelihood
iteration <- iteration + 1
}
# Compute ICs
df <- 2*length(alpha)-1+1+1+1
aic <- -2*loglikelihood + 2*df
bic <- -2*loglikelihood + log(length(lower))*df
return(list(pi = pi, alpha = alpha, beta = beta, shape = shape, theta = theta, gamma = gamma,
loglikelihood = loglikelihood, AIC=aic, BIC=bic))
}
# Fit splicing model using EM algorithm for given M and shapes
# Additional input: eps = numerical tolerance for EM algorithm (stopping criterion),
# beta_tol = numerical tolerence to discard beta's
# maxiter = maximal number of EM iterations
.spliceEM_fit_raw <- function(lower, upper, censored, trunclower = 0, tsplice, truncupper = Inf,
pi, theta, shape, beta, gamma,
eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf) {
n <- length(lower)
M <- length(shape)
# Separate uncensored and censored observations
uncensored <- !censored
# Boolean for having uncensored and censored observations
#no_censoring <- (sum(uncensored) != 0)
#censoring <- (sum(censored) != 0)
# Indices for the five different cases
# Some can be numeric(0)!
# t^l <= x_i=l_i=u_i <= t <T
ind1 <- which(uncensored & (upper <= tsplice))
# t^l < t < x_i=l_i=u_i <=T
ind2 <- which(uncensored & (upper > tsplice))
# t^l <= l_i < u_i <= t <T
ind3 <- which(censored & upper <= tsplice)
# t^l < t <= l_i < u_i <=T
ind4 <- which(censored & lower >= tsplice)
# t^l <= l_i < t < u_i <=T
ind5 <- which(censored & lower < tsplice & upper > tsplice)
lower1 <- lower[ind1]
lower2 <- lower[ind2]
lower3 <- lower[ind3]; upper3 <- upper[ind3]
lower4 <- lower[ind4]; upper4 <- upper[ind4]
lower5 <- lower[ind5]; upper5 <- upper[ind5]
# Truncation probabilities (for ME)
t_probs <- pgamma(tsplice, shape, scale=theta) - pgamma(trunclower, shape, scale=theta)
# MLE for alpha, transform beta to alpha
# alpha_tilde = alpha/(F_ME(t)-F_ME(t^l))
alpha_tilde <- beta / t_probs
alpha <- alpha_tilde / sum(alpha_tilde)
# Compute densities and probabilities for observations from the 5 classes
L <- .spliceEM_densprob(lower1=lower1, lower2=lower2, lower3=lower3, lower4=lower4, lower5=lower5,
upper3=upper3, upper4=upper4, upper5=upper5,
pi=pi, shape=shape, theta=theta, alpha_tilde=alpha_tilde, gamma=gamma,
trunclower=trunclower, tsplice=tsplice, truncupper=truncupper)
x1_dens_nosum <- L$x1_dens_nosum
x1_dens <- L$x1_dens
x2_dens <- L$x2_dens
c3_probs_nosum <- L$c3_probs_nosum
c3_probs <- L$c3_probs
c4_probs <- L$c4_probs
c5_probs_nosum <- L$c5_probs_nosum
c5_probs <- L$c5_probs
# Compute log-likelihood
loglikelihood <- .splice_loglikelihood(x1_dens=x1_dens, x2_dens=x2_dens, c3_probs=c3_probs, c4_probs=c4_probs, c5_probs=c5_probs)
iteration <- 1
old_loglikelihood <- -Inf
while (loglikelihood - old_loglikelihood > eps & iteration <= maxiter) {
old_loglikelihood <- loglikelihood
############
# E-step
############
if (length(ind5) != 0) {
# Compute P(X_i<=t | t^l<=l_i<=t<u_i, Theta^{(h-1)}) and P(X_i>t | t^l<=l_i<=t<u_i, Theta^{(h-1)})
probs <- .spliceEM_probs(lower5=lower5, upper5=upper5, trunclower=trunclower, tsplice=tsplice, truncupper=truncupper,
pi=pi, theta=theta, shape=shape, alpha=alpha, gamma=gamma)
P1 <- probs$P1
P2 <- probs$P2
} else {
P1 <- P2 <- 0
}
##
# pi
n1_h <- length(ind1) + length(ind3) + sum(P1)
n2_h <- n - n1_h
##
# ME
if (length(ind1) != 0) {
# Matrix of z_{ij} with i for rows and j for columns
z1 <- .spliceEM_i_z(x1_dens_nosum=x1_dens_nosum, M=M)
# sum_{i in S_{i.}} x_i
E1_ME <- sum(lower1)
} else {
z1 <- as.matrix(c(0, 0))
E1_ME <- 0
}
if (length(ind3) != 0) {
z3 <- .spliceEM_iii_z(c3_probs_nosum=c3_probs_nosum, M=M)
# sum_{i in S_{iii.}} sum_{j=1}^m z3 * E(X_{ij} | Z_{ij}=1, t^l<=l_i<u_i<=t)
E3_ME <- sum(rowSums(z3 * .spliceEM_Estep_ME_iii(lower3=lower3, upper3=upper3, shape=shape, theta=theta)))
} else {
z3 <- as.matrix(c(0, 0))
E3_ME <- 0
}
if (length(ind5) != 0) {
z5 <- .spliceEM_v_z(c5_probs_nosum=c5_probs_nosum, tsplice=tsplice, alpha_tilde=alpha_tilde,
theta=theta, shape=shape, M=M)
# sum_{i in S_{v.}} P_{1,i} * sum_{j=1}^m z5 * E(X_{ij} | Z_{ij}=1, t^l<=l_i<t<u_i)
E5_ME <- sum(P1 * rowSums(z5 * .spliceEM_Estep_ME_v(lower5=lower5, tsplice=tsplice,
pi=pi, gamma=gamma, shape=shape, theta=theta)))
} else {
z5 <- as.matrix(c(0, 0))
E5_ME <- 0
}
##
# gamma
if (length(ind2) != 0) {
# sum_{i in S_{ii.}} ln(x_i/t)
E2_Pa <- sum(log(lower2/tsplice))
} else {
E2_Pa <- 0
}
if (length(ind4) != 0) {
# sum_{i in S_{iv.}} E( ln(X_i/t) | t^l<t<=l_i<u_i, gamma^{(h-1)})
E4_Pa <- sum(.spliceEM_Estep_Pa_iv(lower4=lower4, upper4=upper4, gamma=gamma, tsplice=tsplice))
} else {
E4_Pa <- 0
}
if (length(ind5) != 0) {
# sum_{i in S_{v.}} P_{2,i} * E( ln(X_i/t) | t^l<l_i<=t<u_i, gamma^{(h-1)})
E5_Pa <- sum(P2 * .spliceEM_Estep_Pa_v(upper5=upper5, gamma=gamma, tsplice=tsplice))
} else {
E5_Pa <- 0
}
############
# M-step
############
##
# pi
pi <- n1_h / n
##
# ME
beta <- (colSums(z1) + colSums(z3) + colSums(z5 * P1)) / n1_h
# Remove small beta's
if (min(beta) < beta_tol) {
shape <- shape[beta > beta_tol]
beta <- beta[beta > beta_tol]
beta <- beta / sum(beta)
}
# Estimate theta
theta <- exp(nlm(.spliceEM_theta_nlm, log(theta), E1_ME=E1_ME, E3_ME=E3_ME, E5_ME=E5_ME,
n1_h=n1_h, beta=beta, shape=shape, trunclower=trunclower, truncupper=tsplice)$estimate)
##
# gamma
gamma_old <- gamma
gamma = 1/n2_h * (E2_Pa + E4_Pa + E5_Pa)
# Solve equation numerically
if (is.finite(truncupper)) {
# Minimise to obtain log(gamma)
.fgamma_trunc <- function(lgammapar, gamma_notrunc) {
# Ensure that gamma is strictly positive
gammapar <- exp(lgammapar)
(gammapar - gamma_notrunc - log(truncupper/tsplice)/((truncupper/tsplice)^(1/gammapar)-1))^2
}
gamma <- exp(nlm(.fgamma_trunc, p=log(gamma_old), gamma_notrunc=gamma)$estimate)
}
#######
iteration <- iteration + 1
# truncation probabilities
t_probs <- pgamma(tsplice, shape, scale=theta) - pgamma(trunclower, shape, scale=theta)
# MLE for alpha, transform beta to alpha
alpha_tilde <- beta / t_probs
alpha <- alpha_tilde / sum(alpha_tilde)
# Compute densities and probabilities for observations from the 5 classes
L <- .spliceEM_densprob(lower1=lower1, lower2=lower2, lower3=lower3, lower4=lower4, lower5=lower5,
upper3=upper3, upper4=upper4, upper5=upper5,
pi=pi, shape=shape, theta=theta, alpha_tilde=alpha_tilde, gamma=gamma,
trunclower=trunclower, tsplice=tsplice, truncupper=truncupper)
x1_dens_nosum <- L$x1_dens_nosum
x1_dens <- L$x1_dens
x2_dens <- L$x2_dens
c3_probs_nosum <- L$c3_probs_nosum
c3_probs <- L$c3_probs
c4_probs <- L$c4_probs
c5_probs_nosum <- L$c5_probs_nosum
c5_probs <- L$c5_probs
# Compute log-likelihood
loglikelihood <- .splice_loglikelihood(x1_dens=x1_dens, x2_dens=x2_dens, c3_probs=c3_probs, c4_probs=c4_probs, c5_probs=c5_probs)
}
# Compute ICs
df <- 2*length(alpha)-1+1+1+1
aic <- -2*loglikelihood + 2*df
bic <- -2*loglikelihood + log(n)*df
return(list(pi = pi, alpha = alpha, beta = beta, shape = shape, theta = theta, gamma = gamma,
loglikelihood = loglikelihood,
iteration = iteration, AIC = aic, BIC = bic))
}
# Compute densities and probabilities for observations from the 5 classes
.spliceEM_densprob <- function(lower1, lower2, lower3, lower4, lower5, upper3, upper4, upper5,
pi, shape, theta, alpha_tilde, gamma,
trunclower, tsplice, truncupper) {
n1 <- length(lower1)
n3 <- length(lower3)
n5 <- length(lower5)
L <- list()
if (n1 != 0) {
# Matrix containing Erlang densities for uncensored observations in ME part and all M values of shape
# Note F_1(x) = sum_{j=1}^M beta_j F_E^t(x) = sum_{j=1}^M beta_j F_E(x)/(F_E(t)-F_E(t^l))
# and alpha_tilde = beta/(F_E(t)-F_E(t^l))
L$x1_dens_nosum <- t(t(outer(lower1, shape, dgamma, scale=theta))*alpha_tilde)
# Combine into density vector (multiply with pi because splicing density)
L$x1_dens <- pi * rowSums(L$x1_dens_nosum)
} else {
L$x1_dens_nosum <- as.matrix(c(0, 0))
L$x1_dens <- numeric(0)
}
# Densities of uncensored observations in EVT part
# No problem when ind2=numeric(0)
L$x2_dens <- (1-pi) * .dtpareto(lower2, gamma=gamma, tsplice=tsplice, truncupper=truncupper)
if (n3 != 0) {
# Matrix containing Erlang probabilities for censored observations of type 3 and all M values of shape
L$c3_probs_nosum <- t(t(outer(upper3, shape, pgamma, scale=theta))*alpha_tilde) - t(t(outer(lower3, shape, pgamma, scale=theta))*alpha_tilde)
# Combine into probability vector (multiply with pi because splicing density): F(u)-F(l)
L$c3_probs <- pi * rowSums(L$c3_probs_nosum)
} else {
# Problem when only right-censoring
L$c3_probs_nosum <- as.matrix(c(0, 0))
L$c3_probs <- numeric(0)
}
# Vector of probabilities for uncensored observations of type 4: F(u)-F(l)
# No problem when ind4=numeric(0)
L$c4_probs <- pi + (1-pi) * (.ptpareto(upper4, gamma=gamma, tsplice=tsplice, truncupper=truncupper)-
.ptpareto(lower4, gamma=gamma, tsplice=tsplice, truncupper=truncupper))
if (n5 != 0) {
# Matrix containing Erlang probabilities for censored observations of type 5 and all M values of shape
L$c5_probs_nosum <- t(t(outer(lower5, shape, pgamma, scale=theta))*alpha_tilde)
# Vector of probabilities for uncensored observations of type 5: F(u)-F(l)
L$c5_probs <- pi + (1-pi) * .ptpareto(upper5, gamma=gamma, tsplice=tsplice, truncupper=truncupper) - pi * rowSums(L$c5_probs_nosum)
} else {
# Problem when no intervals over splicing point
L$c5_probs_nosum <- as.matrix(c(0, 0))
L$c5_probs <- numeric(0)
}
return(L)
}
# Splicing log-likelihood
.splice_loglikelihood <- function(x1_dens, x2_dens, c3_probs, c4_probs, c5_probs) {
likelihood_contribution <- numeric(0)
# likelihood contribution (uncensored)
likelihood_contribution <- c(x1_dens, x2_dens)
# likelihood contribution (censored)
likelihood_contribution <- c(likelihood_contribution, c3_probs, c4_probs, c5_probs)
# Compute log-likelihood and put negative contributions equal to -1000
# ind0 <- which(likelihood_contribution<=0)
# loglikelihood_contribution <- likelihood_contribution
# loglikelihood_contribution[!ind0] <- log(likelihood_contribution[!ind0])
# loglikelihood_contribution[ind0] <- -1000
loglikelihood_contribution <- ifelse(likelihood_contribution > 0, log(likelihood_contribution), -1000)
# log-likelihood
return(sum(loglikelihood_contribution))
}
####################
# E-step auxiliary functions
# Compute P(X_i<=t | t^l<=l_i<=t<u_i, Theta^{(h-1)}) and P(X_i>t | t^l<=l_i<=t<u_i, Theta^{(h-1)})
# Use pi^{(h-1)}, gamma^{(h-1)}, ...
.spliceEM_probs <- function(lower5, upper5, trunclower, tsplice, truncupper, pi, theta, shape, alpha, gamma) {
# CDF of ME evaluated in l_i (for all i with t^l<=l_i<=t<u_i)
p1 <- .pME(lower5, trunclower=trunclower, truncupper=tsplice, theta=theta, shape=shape, alpha=alpha)
P1 <- pi*(1-p1) / (pi+(1-pi)*.ptpareto(upper5, gamma=gamma, tsplice=tsplice, truncupper=truncupper) - pi*p1)
P2 <- 1-P1
return(list(P1=P1, P2=P2))
}
# ^{i.}z_{ij}^{(h)}: posterior probabilities (uncensored)
.spliceEM_i_z <- function(x1_dens_nosum, M) {
z1 <- x1_dens_nosum / rowSums(x1_dens_nosum)
# in case all ^{u}z_{ij}^{(h)} for j=1,...,M are numerically 0
z1[is.nan(z1)] <- 1/M
return(z1)
}
# ^{iii.}z_{ij}^{(h)}: posterior probabilities (censored)
.spliceEM_iii_z <- function(c3_probs_nosum, M) {
z3 <- c3_probs_nosum / rowSums(c3_probs_nosum)
# in case all ^{c}z_{ij}^{(h)} for j=1,...,M are numerically 0
z3[is.nan(z3)] <- 1/M
return(z3)
}
# ^{v.}z_{ij}^{(h)}: posterior probabilities (censored)
.spliceEM_v_z <- function(c5_probs_nosum, tsplice, alpha_tilde, shape, theta, M) {
# Use alpha_tilde since also used in c5_probs_nosum
pt <- alpha_tilde*pgamma(tsplice, shape=shape, scale=theta)
z5 <- -sweep(c5_probs_nosum, 2, pt , "-")
z5 <- z5 / rowSums(z5)
# in case all ^{cc}z_{ij}^{(h)} for j=1,...,M are numerically 0
z5[is.nan(z5)] <- 1/M
return(z5)
}
## Expected value of censored observations with l_i<u_i<t for ME
.spliceEM_Estep_ME_iii <- function(lower3, upper3, shape, theta) {
# E3_ME is a matrix!
E3_ME <- theta * (outer(upper3, shape+1, pgamma, scale=theta) - outer(lower3, shape+1, pgamma, scale=theta)) / (outer(upper3, shape, pgamma, scale=theta) - outer(lower3, shape, pgamma, scale=theta))
E3_ME <- t(t(E3_ME)*shape)
# replace numerical 0/0 (NaN) or Inf by correct expected value
E3_ME <- ifelse(is.nan(E3_ME) | is.infinite(E3_ME), ifelse(outer(lower3, shape*theta, ">"), lower3, upper3), E3_ME)
return(E3_ME)
}
## Expected value of censored observations with l_i<t<u_i for ME
.spliceEM_Estep_ME_v <- function(lower5, tsplice, pi, gamma, shape, theta) {
# ## E5_ME is a matrix!
# E5_ME <- .spliceEM_Estep_ME_c(upper3=rep(tsplice,length(lower5)), lower3=lower5, shape=shape, theta=theta)
# E5_ME <- sweep(E5_ME, 1, (1-pi) * (upper5^(-1/gamma+1)-tsplice^(-1/gamma+1))/((gamma-1)*tsplice^(-1/gamma)), "+")
E5_ME <- .spliceEM_Estep_ME_iii(lower3=lower5, upper3=rep(tsplice, length(lower5)), shape=shape, theta=theta)
return(E5_ME)
}
# E-step for Pareto: part of E( ln(X_i/t) | t^l<t<=l_i<u_i, gamma^{(h-1)})
.spliceEM_Estep_Pa_iv <- function(lower4, upper4, gamma, tsplice) {
# Problem if upper4 is infinite
u4 <- ifelse(is.finite(upper4), upper4/tsplice, 1)
u4_gamma <- ifelse(is.finite(upper4), (upper4/tsplice)^(-1/gamma), 0)
return( ((log(lower4/tsplice)+gamma)*(lower4/tsplice)^(-1/gamma) - (log(u4)+gamma)*u4_gamma) / ((lower4/tsplice)^(-1/gamma) - u4_gamma) )
}
# E-step for Pareto: part of E( ln(X_i/t) | t^l<l_i<=t<u_i, gamma^{(h-1)})
.spliceEM_Estep_Pa_v <- function(upper5, gamma, tsplice) {
return( .spliceEM_Estep_Pa_iv(lower4=tsplice, upper4=upper5, gamma=gamma, tsplice=tsplice) )
}
####################
# M-step auxiliary functions
# Auxiliary functions used to estimate theta in the M-step
# This function returns log of optimal theta!
.spliceEM_theta_nlm <- function(ltheta, E1_ME, E3_ME, E5_ME, n1_h, beta, shape, trunclower, truncupper) {
theta <- exp(ltheta)
TT <- .spliceEM_T(trunclower=trunclower, truncupper=truncupper, shape=shape, theta=theta, beta=beta)
return( (theta - ((E1_ME + E3_ME + E5_ME)/n1_h-TT) / sum(beta*shape)) ^ 2 )
}
.spliceEM_T <- function(trunclower, truncupper, shape, theta, beta) {
# Avoid NaN
if (truncupper == Inf) {
# Take log first for numerical stability (avoid Inf / Inf)
deriv_trunc_log_1 <- shape*log(trunclower)-trunclower/theta - (shape-1)*log(theta) - lgamma(shape) - log(1 - pgamma(trunclower, shape, scale=theta))
deriv_trunc <- exp(deriv_trunc_log_1)
} else {
deriv_trunc_log_1 <- shape*log(trunclower)-trunclower/theta - (shape-1)*log(theta) - lgamma(shape) - log(pgamma(truncupper, shape, scale=theta) - pgamma(trunclower, shape, scale=theta))
deriv_trunc_log_2 <- shape*log(truncupper)-truncupper/theta - (shape-1)*log(theta) - lgamma(shape) - log(pgamma(truncupper, shape, scale=theta) - pgamma(trunclower, shape, scale=theta))
deriv_trunc <- exp(deriv_trunc_log_1)-exp(deriv_trunc_log_2)
}
return(sum(beta*deriv_trunc))
}
#####################
# Fast CDFs and PDFs
# ME CDF (faster version based on .ME_cdf)
.pME <- function(x, theta, shape, alpha, trunclower = 0, truncupper = Inf) {
# Matrix with pgamma(x, shape, scale=theta) for elements of x (rows) and elements of shape (columns)
cdf <- outer(x, shape, pgamma, scale=theta)
# Compute CDF as product with alpha and then sum per column
p <- colSums(t(cdf)*alpha)
if (!(trunclower == 0 & truncupper == Inf)) {
l <- .pME(trunclower, theta=theta, shape=shape, alpha=alpha)
u <- .pME(truncupper, theta=theta, shape=shape, alpha=alpha)
p <- ((p - l) / (u - l)) ^ {(x <= truncupper)} * (trunclower <= x)
}
return(p)
}
# Pareto PDF
.dpareto <- function(x, gamma, tsplice) {
ifelse(x <= tsplice, 0, 1/(gamma*tsplice) * (x/tsplice)^(-1/gamma-1))
}
# Pareto CDF
.ppareto <- function(x, gamma, tsplice) {
ifelse(x <= tsplice, 0, 1-(x/tsplice)^(-1/gamma))
}
# Truncated Pareto PDF
.dtpareto <- function(x, gamma, tsplice, truncupper) {
ifelse(x <= truncupper, .dpareto(x=x, gamma=gamma, tsplice=tsplice)/.ppareto(x=truncupper, gamma=gamma, tsplice=tsplice), 0)
}
# Truncated Pareto CDF
.ptpareto <- function(x, gamma, tsplice, truncupper) {
ifelse(x <= truncupper, .ppareto(x=x, gamma=gamma, tsplice=tsplice)/.ppareto(x=truncupper, gamma=gamma, tsplice=tsplice), 1)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.