廣義線性模型中的Gauss Seidel 迭代算法實現(xiàn)

數(shù)值模擬的算法迭代公式推導

R代碼實現(xiàn)

根據(jù)以上公式，代入迭代步驟，即可實現(xiàn)算法。

##------數(shù)據(jù)模擬--------

library(MASS)
##mvrnorm()

##定義一個產(chǎn)生多元正態(tài)分布的隨機向量協(xié)方差矩陣
Simu_Multi_Norm<-function(x_len, sd = 1, pho = 0.5){
  #初始化協(xié)方差矩陣
  V <- matrix(data = NA, nrow = x_len, ncol = x_len)
 
  #mean及sd分別為隨機向量x的均值和方差
 
  #對協(xié)方差矩陣進行賦值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)
}

##產(chǎn)生模擬數(shù)值自變量X
set.seed(123)
X<-mvrnorm(n = 200, mu = rep(0,10), Simu_Multi_Norm(x_len = 10,sd  = 1, pho = 0.5))

##產(chǎn)生模擬數(shù)值：響應變量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)

##產(chǎn)生model matrix
mydata<-data.frame(X = X, y = y)

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

##包含截矩項的系數(shù)
#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)
}

##初始化系數(shù)
b0<-rep(0,length(b_real))

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

##b1用于記錄更新系數(shù)
b1<-b0

##b.best用于存放歷史最大似然值對應系數(shù)
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
#初始化一個列表用于存放每一次迭代的系數(shù)結(jié)果
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循環(huán)，對每個系數(shù)最優(yōu)化
    
    #線性部分
    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)
    
    #     #若似然值有所增加，則將當前系數(shù)保存
    if(ll.new > ll.old){
      #更新系數(shù)
      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
    ##當達到停止條件時，跳出循環(huán)
    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)
迭代結(jié)果如下：

迭代206步收斂，系數(shù)結(jié)果非常接近R內(nèi)部函數(shù)glm()運行的結(jié)果，甚至稍好于這一結(jié)果。