posterior {aphid}

Posterior decoding.

Description

Calculate the posterior probability of a sequence given a model.

Usage

posterior(x, y, ...)

## S3 method for class 'HMM'
posterior(x, y, logspace = "autodetect", cpp = TRUE, ...)

## S3 method for class 'PHMM'
posterior(x, y, logspace = "autodetect", cpp = TRUE, ...)

Arguments

x

an object of class 'HMM' or 'PHMM'.

y

a vector of mode "character" or "raw" (a "DNAbin" or "AAbin" object) representing a single sequence hypothetically emitted by the model in x.

...

additional arguments to be passed between methods.

logspace

logical indicating whether the emission and transition probabilities of x are logged. If logspace = "autodetect" (default setting), the function will automatically detect if the probabilities are logged, returning an error if inconsistencies are found. Note that choosing the latter option increases the computational overhead; therefore specifying TRUE or FALSE can reduce the running time.

cpp

logical, indicates whether the dynamic programming matrix should be filled using compiled C++ functions (default; many times faster). The R version is primarily retained for bug-fixing and experimentation.

Details

See Durbin et al (1998) chapter 3.2 for details on the calculation and interpretation of posterior state probabilities. Currently no method is available for profile HMMs, but this may be included in a future version if required.

Value

a vector, matrix or array of posterior probabilities.

Author(s)

Shaun Wilkinson

References

Durbin R, Eddy SR, Krogh A, Mitchison G (1998) Biological sequence analysis: probabilistic models of proteins and nucleic acids. Cambridge University Press, Cambridge, United Kingdom.

See Also

forward, backward, Viterbi

Examples

  ## Posterior decoding for standard hidden Markov models
  ## The dishonest casino example from Durbin et al (1998) chapter 3.2
  states <- c("Begin", "Fair", "Loaded")
  residues <- paste(1:6)
  ### Define the transition probability matrix
  A <- matrix(c(0, 0, 0, 0.99, 0.95, 0.1, 0.01, 0.05, 0.9), nrow = 3)
  dimnames(A) <- list(from = states, to = states)
  ### Define the emission probability matrix
  E <- matrix(c(rep(1/6, 6), rep(1/10, 5), 1/2), nrow = 2, byrow = TRUE)
  dimnames(E) <- list(states = states[-1], residues = residues)
  ### Build and plot the HMM object
  x <- structure(list(A = A, E = E), class = "HMM")
  plot(x, main = "Dishonest casino HMM")
  ### Calculate posterior probabilities
  data(casino)
  casino.post <- posterior(x, casino)
  plot(1:300, casino.post[1, ], type = "l", xlab = "Roll number",
       ylab = "Posterior probability of dice being fair",
       main = "The dishonest casino")