epiGWAS (version 1.0.2)

forward_sample: Applies the forward algorithm to a single observation

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

Value

Joint probability for the state x in a log form

Details

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

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

# NOT RUN {
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)

# }