fanshenzhang
您没有写例子,那我就帮您写个最小化例子吧。
n = 150
x1 = rnorm(n, 7, 2)
x2 = rnorm(n, 9, 3)
y1 = rnorm(n, 2*x1 + 3*x2 +1, 1)
model1.lm <- lm(y1 ~ x1 +x2)
summary(model1.lm)
输出结果:
Call:
lm(formula = y1 ~ x1 + x2)
Residuals:
Min 1Q Median 3Q Max
-1.94609 -0.63838 0.03421 0.63044 2.59260
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.10162 0.40176 2.742 0.00687 **
x1 1.98770 0.04632 42.909 < 2e-16 ***
x2 2.99571 0.02528 118.520 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1 on 147 degrees of freedom
Multiple R-squared: 0.9908, Adjusted R-squared: 0.9907
F-statistic: 7928 on 2 and 147 DF, p-value: < 2.2e-16
下面是用极大似然法估计:
library(bbmle)
linear_reg_fun <- function(a, b1, b2, sigma){
Y.pred = a + b1*x1 +b2*x2
return(-sum(dnorm(y1, mean = Y.pred, sd = sigma, log = TRUE)) )
}
mle_model <- bbmle::mle2(linear_reg_fun,
start = list(a = 1, b1 = 2, b2 = 3, sigma = 1))
summary(mle_model)
用 MLE 的输出结果,可以看出两种方法估值非常相近:
Maximum likelihood estimation
Call:
bbmle::mle2(minuslogl = linear_reg_fun, start = list(a = 1, b1 = 2,
b2 = 3, sigma = 1))
Coefficients:
Estimate Std. Error z value Pr(z)
a 1.101617 0.397719 2.7698 0.005608 **
b1 1.987698 0.045858 43.3445 < 2.2e-16 ***
b2 2.995711 0.025022 119.7235 < 2.2e-16 ***
sigma 0.990272 0.057173 17.3205 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-2 log L: 422.7485
上述的例子是线性回归方程,您可以用这个方法套用到估计 非线性多元回归参数的例子里。