数值模拟的算法迭代公式推导
R代码实现
根据以上公式,代入迭代步骤,即可实现算法。
library(MASS)
Simu_Multi_Norm<-function(x_len, sd = 1, pho = 0.5){
V <- matrix(data = NA, nrow = x_len, ncol = x_len)
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)
}
set.seed(123)
X<-mvrnorm(n = 200, mu = rep(0,10), Simu_Multi_Norm(x_len = 10,sd = 1, pho = 0.5))
beta<-c(1,2,0,0,3,0,0,0,-2,0)
prob<-exp( X %*% beta)/(1+exp( X %*% beta))
y<-rbinom(n = 200, size = 1,p = prob)
mydata<-data.frame(X = X, y = y)
b_real<-beta
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))
b1<-b0
b.best<-b0
ll.old<-loglikelihood(X = X,y = y, b = b0)
diff<-1
iter<-0
epsi<-1e-10
max_iter<-10000
b_history<-list(data.frame(b0))
ll_list<-list(ll.old)
while(diff > epsi & iter < max_iter){
for(j in 1:length(b_real)){
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)
ll_hist<-do.call(rbind,ll_list)
iter
ll.best<-max(ll_hist)
ll.best
b.best
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)
迭代结果如下:
迭代206步收敛,系数结果非常接近R内部函数glm()运行的结果,甚至稍好于这一结果。