PELT.R

PELT <- function(data, measure, beta, K){
  n = length(data)
  F = rep(0, n + 1); F[1] = -beta;
  cp = rep(NULL, n)
  R = rep(list(c(0)), n+1)
  for (t in 2:(n+1)){
    hold = rep(0, length(R[[t]]))
    for (j in R[[t]]) {hold = F[j] + measure(data[(j+1):t])}
    F[t] = min(hold)
    tau = R[[t]][which.min(hold)]
    cp[t] = list(cp[tau], tau)
    
    for (tau in R[[t]]){
      if (F(tau) + measure(data[(tau + 1) : t]) + K <= Fsample(base, 1000, replace = TRUE, prob = c(0.1, 0.1, 0.4, 0.4))[t]){
        R[[t + 1]] = c(R[[t + 1]], tau)
      }
    }
    R[[t + 1]] = c(R[[t + 1]], t)
  }
  return (cp[n])
}

measure = loglikhood

base = c("a", "t", "c", "g")
samp = sample(base, 1000, replace = TRUE, prob = c(0.1, 0.1, 0.4, 0.4))
samp = c(samp, sample(base, 1000, replace = TRUE, prob = c(0.4, 0.4, 0.1, 0.1)))

PELT(samp, loglikhood, beta = 0.1, K = 0.1)


R = list(c(0))
c(R,c(1))

rep(list(c(0)),100)

loglikhood(samp)
library(markv)
aaaaana/changepoints documentation built on May 26, 2019, 1:35 p.m.