广义线性模型中的Gauss Seidel 迭代算法实现

浏览: 2041

数值模拟的算法迭代公式推导

Clipboard Image.png

Clipboard Image.png

R代码实现

根据以上公式,代入迭代步骤,即可实现算法。

##------数据模拟--------

library(MASS)
##mvrnorm()

##定义一个产生多元正态分布的随机向量协方差矩阵
Simu_Multi_Norm<-function(x_len, sd = 1, pho = 0.5){
#初始化协方差矩阵
V <- matrix(data = NA, nrow = x_len, ncol = x_len)

#mean及sd分别为随机向量x的均值和方差

#对协方差矩阵进行赋值pho(i,j) = pho^|i-j|
for(i in 1:x_len){ ##遍历每一行
for(j in 1:x_len){ ##遍历每一列
V[i,j] <- pho^abs(i-j)
}
}

V<-(sd^2) * V
return(V)
}


##产生模拟数值自变量X
set.seed(123)
X<-mvrnorm(n = 200, mu = rep(0,10), Simu_Multi_Norm(x_len = 10,sd = 1, pho = 0.5))

##产生模拟数值:响应变量y
beta<-c(1,2,0,0,3,0,0,0,-2,0)

#alpha<-0

#prob<-exp(alpha + X %*% beta)/(1+exp(alpha + X %*% beta))

prob<-exp( X %*% beta)/(1+exp( X %*% beta))

y<-rbinom(n = 200, size = 1,p = prob)

##产生model matrix
mydata<-data.frame(X = X, y = y)

#X<-model.matrix(y~., data = mydata)



##包含截矩项的系数
#b_real<-c(alpha,beta)
b_real<-beta

#define the log-likelihood function
loglikelihood<-function(X, y, b){
linear_comb<-as.vector(X %*% b)
ll<-sum(y*linear_comb) + sum(log(1/(1+exp(linear_comb))))
return (ll)
}


##初始化系数
b0<-rep(0,length(b_real))

#b0<- b_real+rnorm(length(b_real), mean = 0, sd = 0.1)


##b1用于记录更新系数
b1<-b0

##b.best用于存放历史最大似然值对应系数
b.best<-b0

# the initial value of loglikelihood

ll.old<-loglikelihood(X = X,y = y, b = b0)


# initialize the difference between the two steps of theta
diff<-1
#record the number of iterations
iter<-0
#set the threshold to stop iterations
epsi<-1e-10
#the maximum iterations
max_iter<-10000
#初始化一个列表用于存放每一次迭代的系数结果
b_history<-list(data.frame(b0))

#初始化列表用于存放似然值
ll_list<-list(ll.old)



#-------Gauss-Seidel 迭代-------
while(diff > epsi & iter < max_iter){
for(j in 1:length(b_real)){
#对j循环,对每个系数最优化

#线性部分
linear_comb<-as.vector(X %*% b0)

#分子
nominator<-sum(y*X[,j] - X[,j] * exp(linear_comb)/(1+exp(linear_comb)))
#分母,即二阶导部分
denominator<- -sum(X[,j]^2 * exp(linear_comb)/(1+exp(linear_comb))^2)
#
b0[j]<-b0[j] - nominator/denominator
#更新似然值
ll.new<- loglikelihood(X = X, y = y, b = b0)

# #若似然值有所增加,则将当前系数保存
if(ll.new > ll.old){
#更新系数
b.best[j]<-b0[j]
}

#求差异
diff<- abs((ll.new - ll.old)/ll.old)
ll.old <- ll.new
iter<- iter+1
b_history[[iter]]<-data.frame(b0)
ll_list[[iter]]<-ll.old
##当达到停止条件时,跳出循环
if(diff < epsi){
break
}

}


}

b_hist<-do.call(rbind,b_history)
#b_hist

ll_hist<-do.call(rbind,ll_list)
#ll_hist

#
iter

##
ll.best<-max(ll_hist)
ll.best

##
b.best




##---------glm()验证-------

my_glm<-glm(y~0 + X.1 + X.2 + X.3+ X.4+ X.5+ X.6+ X.7+ X.8+ X.9+ X.10,
data = mydata,family = binomial(link = "logit"))
summary(my_glm)

coeff_glm<-my_glm$coefficients

cbind(coeff_glm,b.best,b_real)

迭代结果如下:


Clipboard Image.png

迭代206步收敛,系数结果非常接近R内部函数glm()运行的结果,甚至稍好于这一结果。

推荐 0
本文由 黄耀鹏 创作,采用 知识共享署名-相同方式共享 3.0 中国大陆许可协议 进行许可。
转载、引用前需联系作者,并署名作者且注明文章出处。
本站文章版权归原作者及原出处所有 。内容为作者个人观点, 并不代表本站赞同其观点和对其真实性负责。本站是一个个人学习交流的平台,并不用于任何商业目的,如果有任何问题,请及时联系我们,我们将根据著作权人的要求,立即更正或者删除有关内容。本站拥有对此声明的最终解释权。

0 个评论

要回复文章请先登录注册