2n_ril.R

# 2^n-way RIL by selfing
# binned two-locus probabilities, likelihood for rec frac, MLE

# 8-way RIL by selfing

probs8 <- function(r) c((1-r)^2, r*(1-r), r/2)/8/(1+2*r) * c(8, 8, 48)


llik8 <-
    function(r, counts)
{
    u <- counts[1]
    v <- counts[2]
    w <- counts[3]
    n <- sum(counts)

    2*u*log(1-r) + v*log(r*(1-r)) + w*log(r) - n*log(1+2*r)
}

mle8 <- function(counts)
{
    u <- counts[1]
    v <- counts[2]
    w <- counts[3]
    n <- sum(counts)

    A <- sqrt(4*n^2 + 4*n*(2*u - 2*v - 3*w) + 9*w^2 + 12*w*(u+2*v) + 16*v^2 + 16*u*v + 4*u^2)

    (2*n + 2*u - w - A)/4/(n - w - 2*v - 2*u)
}


# 16-way RIL by selfing
probs16 <- function(r) c((1-r)^3, r*(1-r)^2, r*(1-r), r)/c(16,16,32,64)/(1+2*r) * c(16, 16, 32, 192)

llik16 <-
    function(r, counts)
{
    u <- counts[1]
    v <- counts[2]
    w <- counts[3]
    y <- counts[4]
    n <- sum(counts)

    (3*u+2*v)*log(1-r) + (v+w+y)*log(r) - n*log(1+2*r)
}

mle16 <- function(counts)
{
    u <- counts[1]
    v <- counts[2]
    w <- counts[3]
    y <- counts[4]
    n <- sum(counts)

    A <- sqrt(9*y^2 + 6*y*(3*w + 5*v + 3*u - 2*n) + 9*w^2 + 6*w*(5*v+3*u-2*n) + 25*v^2 + 2*v*(15*u-2*n) + 9*u^2 + 12*n*u +4*n^2)

    (A + y + w - v - 3*u - 2*n)/4/(y + w + 3*v + 3*u - n)

}