在R中编写正确正态对数似然函数及参数推断问题求助

我来帮你梳理并修正代码里的问题，你的原函数在似然表达式、参数处理上都有小疏漏，咱们逐个搞定：

首先：原代码的核心问题 #

1. 似然公式写错了：你的模型里每个观测的方差是tau² + u[i]²（因为u是已知向量，每个样本对应不同的u值），所以不能用n*log(tau²+u²)，得对每个i计算log(tau² + u[i]²)再求和；另外平方项的求和应该是sum( (x - mu)² / (tau² + u²) )，原代码里的2*tau²+u²是错误的（正态对数似然里是1/(2*方差)，你把系数和方差混在一起了）。
2. sapply用错了：对数似然函数对一组参数（mu和tau各一个值）应该返回一个标量，不需要用sapply遍历c(mu,tau)，这会导致函数返回向量而非单个值，后续optim也无法正常工作。
3. 参数传递失误：sapply里的匿名函数没有正确接收mu和tau参数，导致计算时变量未定义。

第一步：编写正确的对数似然函数 #

先把已知的u向量定义好，然后写出符合模型的对数似然：

# 给定的数据
x <- c(3.3569,1.9247,3.6156,1.8446,2.2196,6.8194,2.0820,4.1293,0.3609,2.6197)
u <- c(1,3,0.5,0.2,2,1.7,0.4,1.2,1.1,0.7)

# 正确的对数似然函数
normal.lik <- function(theta, x, u) {
  mu <- theta[1]
  tau <- theta[2]
  
  # 每个观测对应的方差
  variances <- tau^2 + u^2
  n <- length(x)
  
  # 计算对数似然
  logl <- -0.5 * n * log(2 * pi) - 0.5 * sum(log(variances)) - 0.5 * sum( (x - mu)^2 / variances )
  
  return(logl)
}

测试一下这个函数，传入参数c(1,2)：

normal.lik(c(1,2), x, u)
# 会返回一个标量，比如我的运行结果是 -28.70475

第二步：固定tau，绘制mu的对数似然曲线 #

你已经定义了mu的序列，现在固定tau=2，计算每个mu对应的对数似然并绘图：

mu <- seq(0,10,length=1000)
tau_fixed <- 2

# 计算每个mu对应的对数似然值
logl_vals <- sapply(mu, function(m) normal.lik(c(m, tau_fixed), x, u))

# 绘制曲线
plot(mu, logl_vals, type = "l", lwd = 2, xlab = expression(mu), ylab = "对数似然值", main = "固定tau=2时，mu的对数似然曲线")
abline(v = mu[which.max(logl_vals)], col = "red", lty = 2) # 标记似然最大值对应的mu

运行这段代码就能得到你想要的曲线，红色虚线就是tau=2时mu的MLE。

第三步：用optim求解mu和tau的MLE #

optim默认求最小值，所以我们可以通过control=list(fnscale=-1)让它最大化对数似然。另外，tau是方差的组成部分，必须大于0，所以可以加个约束（比如把tau参数化为指数形式，保证非负），这样结果更可靠：

方法1：直接最大化对数似然 #

# 初始参数值
start_theta <- c(1, 1)

# 求解MLE
theta_hat <- optim(start_theta, normal.lik, control = list(fnscale = -1), x = x, u = u)$par
names(theta_hat) <- c("mu的MLE", "tau的MLE")

# 输出结果
theta_hat

方法2：约束tau>0（更严谨） #

把tau用指数形式表示，确保它始终为正：

normal.lik_constrained <- function(theta, x, u) {
  mu <- theta[1]
  tau <- exp(theta[2]) # 用指数转换保证tau>0
  
  variances <- tau^2 + u^2
  n <- length(x)
  
  logl <- -0.5 * n * log(2 * pi) - 0.5 * sum(log(variances)) - 0.5 * sum( (x - mu)^2 / variances )
  
  return(logl)
}

# 初始值：theta[2]对应log(tau)，初始tau=1所以theta[2]=0
start_theta_constrained <- c(1, 0)

theta_hat_constrained <- optim(start_theta_constrained, normal.lik_constrained, control = list(fnscale = -1), x = x, u = u)$par
# 转换回tau的原始值
mu_MLE <- theta_hat_constrained[1]
tau_MLE <- exp(theta_hat_constrained[2])

cat("带约束的MLE结果：\nmu =", round(mu_MLE, 4), "\ntau =", round(tau_MLE, 4), "\n")

运行后就能得到mu和tau的最大似然估计值，你可以代入对数似然函数验证是否为最大值点。

内容的提问来源于stack exchange，提问作者Hans Christensen