HMM

#states: Gizli Markov Modeli'nde (HMM) iki durum olduğunu belirten bir vektör (alpha ve beta).
states <- c("alpha", "beta")

#observations: Gözlemlenen dizi (GGCT).
observations <- strsplit("GGCT", "")[[1]]

#start_prob: Her bir durum için başlangıç olasılıklarını logaritmik formda belirten vektör. Burada alpha ve beta durumlarının başlangıç olasılıkları eşit olarak log(1) (0) 
start_prob <- c(alpha = log(1), beta = log(1))  

#trans_prob: Durumlar arası geçiş olasılıklarını logaritmik formda belirten matris. alpha durumundan diğer durumlara geçiş olasılıkları, beta durumundan diğer durumlara geçiş olasılıklarıdır.
trans_prob <- matrix(c(log(9/10), log(1/10),
                       log(1/10), log(9/10)), nrow=2, byrow=TRUE,
                     dimnames=list(states, states))

#emission_prob: Her bir durum için her bir gözlemin yayılma olasılıklarını logaritmik formda belirten matris
emission_prob <- matrix(c(log(2/5), log(2/5), log(1/10), log(1/10),  # alpha
                          log(1/5), log(1/5), log(3/10), log(3/10)),  # beta
                        nrow=2, byrow=TRUE,
                        dimnames=list(states, c("A", "G", "C", "T")))

# Viterbi algoritmasının uygulanması
viterbi <- function(observations, states, start_prob, trans_prob, emission_prob) {
  n <- length(observations)
  m <- length(states)
  delta <- matrix(NA, nrow=m, ncol=n)
  psi <- matrix(NA, nrow=m, ncol=n)
#n: Gözlem dizisinin uzunluğu.
#m: Durum sayısı.
#delta: Maksimum olasılıkların logaritmalarını tutan matris.
#psi: Backtracking matrisidir ve maksimum olasılıkları sağlayan önceki durumları tutar.
  
  # Başlangıç durumu
  for (i in 1:m) {
    delta[i, 1] <- start_prob[states[i]] + emission_prob[states[i], observations[1]]
    psi[i, 1] <- 0
  }
#İlk gözlem için her bir durumun başlangıç olasılığı ve yayılma olasılığı logaritmalarının toplamı hesaplanır ve delta matrisine atanır.
#psi matrisinin ilk sütunu başlangıç durumu olduğundan 0 olarak atanır.
  
  # Recursion
  for (t in 2:n) {
    for (j in 1:m) {
      max_val <- -Inf
      max_state <- NA
      for (i in 1:m) {
        val <- delta[i, t-1] + trans_prob[states[i], states[j]]
        if (val > max_val) {
          max_val <- val
          max_state <- i
        }
      }
      delta[j, t] <- max_val + emission_prob[states[j], observations[t]]
      psi[j, t] <- max_state
    }
  }
#Her zaman adımı (t=2'den n'e kadar) ve her durumu (j) için:
#Her bir önceki durumdan (i) geçiş olasılıkları ve yayılma olasılıkları hesaplanarak maksimum değer bulunur (max_val) ve ilgili durum (max_state) kaydedilir.
#Bu değerler delta ve psi matrislerine atanır.
  
  # Termination
  p_star <- max(delta[, n])
  q_star <- numeric(n)
  q_star[n] <- which.max(delta[, n])
#delta matrisinin son sütununda en yüksek log-olasılık değeri (p_star) ve bu değeri sağlayan durum (q_star[n]) bulunur.
  
  # Backtracking
  for (t in (n-1):1) {
    q_star[t] <- psi[q_star[t+1], t+1]
  }
#Son durumdan başlayarak (q_star[n]), psi matrisini kullanarak önceki durumlar belirlenir ve en olası durum dizisi (state_sequence) oluşturulur. 
  
  state_sequence <- states[q_star]
  return(list(state_sequence=state_sequence, log_prob=p_star))
}

# HMM ile en olası durum dizisini bulma ve olasılığını hesaplama
result <- viterbi(observations, states, start_prob, trans_prob, emission_prob)

# Sonuçları ekrana yazdırma
cat("En olası durum dizisi:", paste(result$state_sequence, collapse=" "), "\n")

## En olası durum dizisi: beta beta beta beta

cat("En yüksek log-olasılık:", result$log_prob, "\n")

## En yüksek log-olasılık: -5.942903

library(HMM)

## 
## Attaching package: 'HMM'

## The following object is masked _by_ '.GlobalEnv':
## 
##     viterbi

# Initialize the HMM
hmm <- initHMM(
  States = c("alpha", "beta"),
  Symbols = c("A", "T", "G", "C"),
  transProbs = matrix(c(log(9/10), log(1/10), log(1/10), log(9/10)), 2),
  emissionProbs = matrix(c(
    log(0.4), log(0.1), log(0.4), log(0.1),
    log(0.2), log(0.3), log(0.2), log(0.3)
  ), nrow = 2, byrow = TRUE)
)

print(hmm)

## $States
## [1] "alpha" "beta" 
## 
## $Symbols
## [1] "A" "T" "G" "C"
## 
## $startProbs
## alpha  beta 
##   0.5   0.5 
## 
## $transProbs
##        to
## from         alpha       beta
##   alpha -0.1053605 -2.3025851
##   beta  -2.3025851 -0.1053605
## 
## $emissionProbs
##        symbols
## states           A         T          G         C
##   alpha -0.9162907 -2.302585 -0.9162907 -2.302585
##   beta  -1.6094379 -1.203973 -1.6094379 -1.203973

start_prob <- c(log(0.5), log(0.5))

 
viterbi <- function(hmm, observations, start_prob) {
  n_states <- length(hmm$States)
  n_obs <- length(observations)
  
  # Initialization
  V <- matrix(NA, nrow = n_states, ncol = n_obs)
  path <- matrix(NA, nrow = n_states, ncol = n_obs)
  
  for (s in 1:n_states) {
    V[s, 1] <- start_prob[s] + hmm$emissionProbs[s, observations[1]]
    path[s, 1] <- s
  }
  
  # Recursion
  for (t in 2:n_obs) {
    for (s in 1:n_states) {
      max_prob <- -Inf
      best_state <- NA
      for (s_prev in 1:n_states) {
        prob <- V[s_prev, t-1] + hmm$transProbs[s_prev, s] + hmm$emissionProbs[s, observations[t]]
        if (prob > max_prob) {
          max_prob <- prob
          best_state <- s_prev
        }
      }
      V[s, t] <- max_prob
      path[s, t] <- best_state
    }
  }
  
  # Termination
  best_path <- rep(NA, n_obs)
  best_prob <- -Inf
  for (s in 1:n_states) {
    if (V[s, n_obs] > best_prob) {
      best_prob <- V[s, n_obs]
      best_path[n_obs] <- s
    }
  }
  
  for (t in (n_obs-1):1) {
    best_path[t] <- path[best_path[t+1], t+1]
  }
  
  best_path <- hmm$States[best_path]
  return(list(path = best_path, prob = best_prob))
}

observations <- c("G", "G", "C", "T")
viterbi_result <- viterbi(hmm, observations, start_prob)
print(viterbi_result)

## $path
## [1] "beta" "beta" "beta" "beta"
## 
## $prob
## [1] -6.63605