R/MLEModels.R

In radamsRHA/Misc.Statistics.Genetics.Models:

####################################################################
################# Transitional Probality Matrix; t >> ##############
####################################################################

JC69.Pt<- function(sub.rate, t){
  instant.rate <- sub.rate/3
  p0 <- (1/4)+(3/4)*exp(-4*instant.rate*t)
  p1 <- (1/4)-(1/4)*exp(-4*instant.rate*t)
  Pt <- matrix(nrow = 4, ncol = 4)
  Pt[,1] <- c(p0, p1, p1, p1)
  Pt[,2] <- c(p1, p0, p1, p1)
  Pt[,3] <- c(p1, p1, p0, p1)
  Pt[,4] <- c(p1, p1, p1, p0)
  colnames(Pt) <- c("T", "C", "A", "G")
  rownames(Pt) <- c("T", "C", "A", "G")
  return(Pt)
}

####################################################################
################# Untransformed ####################################
####################################################################

MLE.d.JC69 <- function(d, x, n){
  p1 <- (1/4)-(1/4)*exp(-4*d/3)
  p0 <- 1-3*p1
  MLE.d <- ((1/4)*(p1)^x)*(((1/4)*p0)^(n-x))
  return(MLE.d)
}

####################################################################
################# Log transformed ##################################
####################################################################


MLE.d.JC69 <- function(d, x, n){
  p1 <-(1/4)-(1/4)*exp(-4*d/3)
  p0 <- 1-3*p1
  MLE.d <- x*log((1/4)*p1)+(n-x)*log((1/4)*p0)
  return(MLE.d)
}

# d will be maximized by d.hat
estimate.d <- function(x, n){
  MLE.d <- (-3/4)*log(1-(4/3)*(x/n))
  return(MLE.d)
}