杂交模型中遗传参数的最大似然估计

from wu2006 Example 2.4

一个位点有2个allele——A和a，在F2子代中可构成3中基因型：AA，Aa和aa，分离比一般为p1:p2:(1-p1-p2)=1:2:1。因此观测值是一个混合分布：

残差服从正态分布。因此：
公式1
其对数似然函数为：

这里
最大似然估计是对数似然函数对参数求导，得：

为了简便有：

公式2

于是有：

采用迭代可以得到估计值。
分别对均值和方差求导有：

将式子等零可得：
公式3
公式4

这2个式子是程序中需要的，也就是说上面的一些推到步骤在编写程序时用不到。
这2个式子比较吓人，不容易看懂，和程序对照着看会更容易懂。下面上程序：

y<-c(79,82,85,87,100,101,102,103,124,125,126,127)
# mean(y)
# sd(y)^2
# var(y-mean(y))
# sum((y-mean(y))^2/(length(y)-1))

n <- length(y)
# 迭代次数
nit <- 20

# 均值初始值
MAA <- matrix(max(y), nit,1)
MAa <- matrix(mean(y), nit,1)
Maa <- matrix(min(y), nit,1)
# 标准差和方差初始值
S <- matrix(sd(y), nit,1)
S2 <- matrix(NA, nit,1)
S2[1] <- 5

# 迭代 注意是从第二次开始的 因为要把这一次的结果赋值给第一次 逐步迭代
for (i in 2:nit) {
  # 公式1
  PAA_temp <- 0.25*dnorm(y, mean = MAA[i-1], sd=S[i-1])
  PAa_temp <- 0.5*dnorm(y, mean = MAa[i-1], sd=S[i-1])
  Paa_temp <- 0.25*dnorm(y, mean = Maa[i-1], sd=S[i-1])
  # 公式2
  PAA <- PAA_temp/(PAA_temp+PAa_temp+Paa_temp)
  PAa <- PAa_temp/(PAA_temp+PAa_temp+Paa_temp)
  Paa <- 1-PAA-PAa
  # 公式3 均值估计值
  MAA[i] <- sum(y*PAA)/sum(PAA)
  MAa[i] <- sum(y*PAa)/sum(PAa)
  Maa[i] <- sum(y*Paa)/sum(Paa)
  # 公式4 方差估计值
  S2[i] <- (sum(PAA*(y-MAA[i])^2)+sum(PAa*(y-MAa[i])^2)+sum(Paa*(y-Maa[i])^2))/(n*sum(PAA+PAa+Paa))
  S[i] <- sqrt(S2[i])
}

# 作图
par(mfrow=c(2,2))
par(mar=c(4,4,1,1))
plot(MAA, type = "l", lwd=2)
plot(MAa, type = "l", lwd=2)
plot(Maa, type = "l", lwd=2)
plot(S2, type = "l", lwd=2)

-end-