R/hmm.vitFUN.rils.R

In MPR.genotyping: Maximum Parsimony of Recombination to Infer Parental Genotypes

hmm.vitFUN.rils <-
function (geno, position, geno.probability, transitionFUN = phy2get.haldane.rils, 
            emissionFUN = makeEmissionFUN(errorRate = 0.01), ...) 
  {
    n.obs <- length(geno)
    n.state <- length(geno.probability)
    psi <- delta <- matrix(0, nrow = n.state, ncol = n.obs)
    n.con <- geno.cr <- numeric(n.obs)
    geno.dis <- abs(diff(position))
    n.con[1] <- 1
    g <- geno[1]
    for (i in 2:n.obs) {
      n.con[i] <- ifelse(geno[i] == g, n.con[i - 1] + 1, 1)
      g <- geno[i]
    }
    for (i in 1:n.state) delta[i, 1] <- log(geno.probability[i]) + 
      emissionFUN(i, geno[1], n.con[1])
    preProb <- numeric(n.state)
    for (t in 2:n.obs) {
      for (j in 1:n.state) {
        for (i in 1:n.state) preProb[i] <- delta[i, t - 1] + 
            log(transitionFUN(i, j, geno.dis[t - 1]))
        psi[j, t] <- which.max(preProb)
        delta[j, t] <- max(preProb) + emissionFUN(j, geno[t], 
                                                  n.con[t])
      }
    }
    geno.cr[n.obs] <- which.max(delta[, n.obs])
    for (t in seq(n.obs - 1, 1, by = -1)) geno.cr[t] <- psi[geno.cr[t + 
                                                                      1], t + 1]
    geno.cr
  }