forward_sample: Applies the forward algorithm to a single observation

In EpiSlim/epiGWAS: Robust Methods for Epistasis Detection

Description

The forward algorithm is applied in order to compute the joint probability for the observation x. For hidden Markov models, the forward algorithm is an attractive option because of its linear complexity in the number of hidden states. However, the complexity becomes quadratic in terms of the dimensionality of the latent space.

Usage

forward_sample(x, p_init, p_trans, p_emit)

Arguments

x	one-sample genotype
p_init	marginal distributions for the first hidden state
p_trans	3D dimensional array for the transition probabilities
p_emit	3D dimensional array for the emission probabilities

Details

Our implementation of the forward algorithm makes use of the LogSumExp transformation for increased numerical stability.

Value

Joint probability for the state x in a log form

References

Rabiner, Lawrence R. 'A tutorial on hidden Markov models and selected applications in speech recognition.' Proceedings of the IEEE 77.2 (1989): 257-286.

Examples

p <- 3 # Number of states
K <- 2 # Dimensionality of the latent space

p_init <- rep(1 / K, K)
p_trans <- array(runif((p - 1) * K * K), c(p - 1, K, K))
# Normalizing the transition probabilities
for (j in seq_len(p - 1)) {
  p_trans[j, , ] <- p_trans[j, , ] / (matrix(rowSums(p_trans[j, , ]), ncol = 1) %*% rep(1, K))
}

p_emit <- array(stats::runif(p * 3 * K), c(p, 3, K))
# Normalizing the emission probabilities
for (j in seq_len(p)) {
  p_emit[j, , ] <- p_emit[j, , ] / (matrix(rep(1, 3), ncol = 1) %*% colSums(p_emit[j, , ]))
}

X <- (runif(p, min = 0, max = 1) < 0.5) + (runif(p, min = 0, max = 1) < 0.5)

# Computing the joint log-probabilities
log_prob <- forward_sample(X, p_init, p_trans, p_emit)

EpiSlim/epiGWAS documentation built on Nov. 19, 2019, 7:15 p.m.