在R中编码高斯对数似然编码

我试图通过编码高斯对数似然来解码optim()来学习R，但在经过数小时的汗水之后，我依然没有达到标准。（这是自学，没有家庭作业。）在R中编码高斯对数似然编码

我下面在写这样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

来源

2017-07-08 invictus

您应该为X添加一个拦截列。您能告诉我们您的起始值是什么吗？（不知道当你有三列X加上一个标准偏差时你如何以5个参数结束......） –

你的方法的基础是完善的，但一些细节是错误的。首先，根据高斯线性模型构造数据是有意义的;例如

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()函数参数化在SD，而不是方差的条款。

然后

> 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估计值不同。

我希望这可以帮助

来源

2017-07-08 02:30:59

在R中编码高斯对数似然编码

回答

相关问题