使用R编写高斯对数似然函数

Question

我正在尝试通过编写高斯对数似然函数来学习R语言，以便使用optim()进行求解。但是，经过几个小时的努力，我仍然偏离了正确的方向。（这不是作业而是自学。）

我按照许多用户编写的函数的惯例编写了一个函数，如loglik <- function(theta, y, x)，其中theta是权重向量beta和sigma的组合，y是输出结果，x是数据。

下面是我的完整代码和模拟数据。运行它，您会发现与lm()相比，我的函数偏离了正确的方向。有人能给我一些关于我错在哪里的想法吗？

# random data
set.seed(111)
y = sample(1:100,100)
x1 = sample(1:100,100)*rnorm(1,0)
x2 = sample(x1)*rnorm(1,0)
x3 = sample(x2)*rnorm(1,0)
dat = data.frame(x1,x2,x3)

# define gaussian log-likelihood
logLik <- function(theta, Y, X){
  X           <- as.matrix(X) # convert data to matrix
  k           <- ncol(X) # get the number of columns (independent vars)
  beta        <- theta[1:k] # vector of weights intialized with starting values
  expected_y  <- X %*% beta  # X is dimension (n x k) and beta is dimension (k x 1)
  sigma2      <- theta[k+1] # should be pulled from the last of the starting values vector
  LL          <- sum(dnorm(Y, mean = expected_y, sd = sigma2, log = T)) # sum of the PDF over each observation
  return(-LL)
}

以下是输出内容：

> optim(logLik, par=starting_values, method="Nelder-Mead", Y=y, X=dat, hessian = T)$par
[1]   0.4832514  -0.2276684  -0.3860800  32.7168490 -38.9030319
> coefficients(lm(y~x1+x2+x3))
(Intercept)          x1          x2          x3 
58.17347451 -0.06587320  0.13001865 -0.03624233 

- invictus

---

你应该在X中添加一个截距列。你能展示一下你的起始值是什么吗？（不确定当你的X有三列加上一个标准偏差时，你是如何得到5个参数的...） - Ben Bolker

Answer 1

你的方法基本正确，但有些细节是错误的。首先，根据高斯线性模型构建数据是有意义的；例如：

set.seed(111)
X <- cbind(1, matrix(rnorm(100*3), 100, 3))
y <- X %*% rep(1, 4) + rnorm(100, 0, 2)

starting.values <- c(1, 1, 1, 1, 2) # actual parameters

# define gaussian log-likelihood
logLik <- function(theta, y, X){
  k           <- ncol(X) # get the number of columns (independent vars)
  beta        <- theta[1:k] # vector of weights intialized with starting values
  expected_y  <- X %*% beta  # X is dimension (n x k) and beta is dimension (k x 1)
  sigma      <- theta[k+1] # should be pulled from the last of the starting values vector
  LL          <- sum(dnorm(y, mean = expected_y, sd = sigma, log = TRUE)) # sum of the PDF over each observation
  return(-LL)
}

顺便提一下，*norm() 函数是以标准差为参数化的，而不是方差。

然后

> optim(logLik, par=starting.values, y=y, X=X, method="BFGS")$par
[1] 1.0471420 1.1411523 0.8167656 0.9840397 1.8910201
Warning message:
In dnorm(y, mean = expected_y, sd = sigma, log = TRUE) : NaNs produced

> summary(lm(y ~ X - 1))

Call:
lm(formula = y ~ X - 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5062 -1.3293  0.1371  1.2057  5.8116 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
X1   1.0471     0.1952   5.365 5.61e-07 ***
X2   1.1412     0.1818   6.275 1.00e-08 ***
X3   0.8168     0.1907   4.282 4.39e-05 ***
X4   0.9840     0.2122   4.638 1.11e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.93 on 96 degrees of freedom
Multiple R-squared:  0.5333,    Adjusted R-squared:  0.5138 
F-statistic: 27.42 on 4 and 96 DF,  p-value: 3.468e-15

注意，method="BFGS" 会发出警告但能产生正确的答案；method="Nelder-Mead" 稍微不太准确。另外，注意通常误差标准差的估算与 ML 估算的不同。

希望这有所帮助。

- John Fox