之前总结了几篇混合线性模型的笔记,继续上传。
这部分主要是介绍如何写出似然函数,通过正态分布,线性回归为例子,并通过R语言编程实现。希望大家可以有所收获。
如何写出似然函数,如何使用R语言编程实现:
当n个观测值相互独立,他们的似然函数(等价于联合密度函数)为:
对似然函数,两边求自然对数:
进一步简化:
normal.lik1<-function(theta,y){
mu<-theta[1]
sigma2<-theta[2]
n<-length(y)
logl<- -.5*n*log(2*pi) -.5*n*log(sigma2) - (1/(2*sigma2))*sum((y-mu)**2) return(-logl)
}
使用R的optim
进行迭代,它的迭代方法有:
mu = 10sigma = 5
y = rnorm(100,mu,va) # n是100,个数
optim
进行迭代求解optim(c(2,1),normal.lik1,y=y)
可以看出,平均数为9.7,方差为21.8
# 将常数项去掉n
ormal.lik2<-function(theta,y){
mu<-theta[1]
sigma2<-theta[2]
n = length(y)
z = (y-mu)/sigma # logl = -0.5*n*log(sigma^2) - (1/(2*sigma^2)*sum((y-mu)^2))
logl = -.5*n*log(sigma2) - (1/(2*sigma2))*sum((y-mu)**2) return(-logl)
}
optim(c(2,1),normal.lik2,y=y)
结果和上面结果一致。
参考
http://www.solinx.co/archives/800
y = ax + b
因此函数为:
set.seed(123)
data.x = rnorm(n = 100)
b0 = 10b1 = 5true.y = b0 + data.x*b1
noise = rnorm(100,0,2)
data.y = true.y + noise
dat = data.frame(mu=1,x=data.x,y=data.y)
head(dat)
lm.logl = function(theta,y,x){
n = length(y)
b0 = theta[1]
b1 = theta[2]
sigma = theta[3]
e = y - b0 - b1*x
logl<- -0.5*n*log(2*pi)- n*log(sigma)- sum(e^2)/(2*sigma^2) return(-logl)
}
optim(c(1,1,1),lm.logl,y=dat$y,x=dat$x)
可以看出,b0的估计值为9.7,b1的估计值为4.8(假定的b0=10,b1=5)。
summary(lm(y~x,data=dat))
Call:
lm(formula = y ~ x, data = dat)
Residuals:
Min 1Q Median 3Q Max
-3.815 -1.367 -0.175 1.161 6.581
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.7944 0.1951 50.2 <2e-16 ***
x 4.8951 0.2138 22.9 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.941 on 98 degrees of freedom
Multiple R-squared: 0.8425, Adjusted R-squared: 0.8409
F-statistic: 524.4 on 1 and 98 DF, p-value: < 2.2e-16
可以看到,评估的b0=9.7, b1 = 4.8,结果一致
对上面的似然函数求偏导
得到的结果和最小二乘法结果一致:
一、确定问题的随机变量类型是离散随机变量还是连续随机变量
二、得出问题的概率分布
三、概率函数转为似然函数
四、似然函数取对数
五、求关于某变量的偏导数
六、解似然方程
https://medium.com/quick-code/maximum-likelihood-estimation-for-regression-65f9c99f815d
https://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/06/lecture-06.pdf
http://www.solinx.co/archives/800
https://www.ime.unicamp.br/~cnaber/optim_1.pdf