MALLOWS模型平均-R代码

（本文仅做学习记录所用。）

我要用我最后的力气喊出那句：“Hansen大佬！！好人一生平安！！”

文章是Hansen大佬2017年的《NOTES AND COMMENTS LEAST SQUARES MODEL A VERAGING》，代码在大佬个人网页。

4.20更新

每一行后面的注释是我写的，希望没写错，为后面的逻辑回归模型平均做铺垫。

1.  
   ## This R function computes the Mallows Model Average (MMA) and
2.  
   ## the Jackknife Model Average (JMA) least-squares estimates.
3.  
   ##
4.  
   ## written by
5.  
   ## Bruce E. Hansen
6.  
   ## Department of Economics
7.  
   ## University of Wisconsin
8.  
   ##
9.  
   ## Format:
10.  
    ##
11.  
    ## result <- gma(y,x,method,subset)
12.  
    ## result$betahat
13.  
    ##
14.  
    ## Inputs:
15.  
    ## y nx1 dependent variable
16.  
    ## x nxp regressor matrix
17.  
    ## method 1x1 set to 1 for Mallows model average estimates
18.  
    ## set to 2 for Jackknife model average estimates
19.  
    ## subset 1x1 set to 1 for pure nested subsets
20.  
    ## set to 2 for all combinations of subsets
21.  
    ## mxp input the (mxp) selection matrix, where m is
22.  
    ## the number of models
23.  
    ## Example:
24.  
    ## Suppose there are 3 candidate models.
25.  
    ## Model 1: y=beta1*x1 beta2*x2 e
26.  
    ## Model 2: y=beta1*x1 beta3*x3 e
27.  
    ## Model 3: y=beta1*x1 beta2*x2 beta4*x4 e
28.  
    ## Then subset <- matrix(c(1,1,1,1,0,1,0,1,0,0,0,1),3,4)
29.  
    ##
30.  
    ## Outputs:
31.  
    ## betahat px1 parameter estimate
32.  
    ## w mx1 weight vector
33.  
    ## yhat nx1 fitted values
34.  
    ## ehat nx1 fitted residuals
35.  
    ## r2 1x1 R-squared
36.  
    ## cn 1x1 Value of Mallows criterion or Cross-Validation criterion
37.  
    ##
38.  
    ## Note:
39.  
    ## Package "quadprog" is required.
40.  
    ## For pure nested subsets, the regressors columns should be ordered, with the
41.  
    ## intercept first and then in order of relevance.
42.  
    ## For all combinations of subsets, p is less than about 20.
43.  
     
44.  
    library(quadprog)
45.  
    gmaN <- function(y,x,method,subset){
46.  
     
47.  
    y <- as.matrix(y)
48.  
    x <- as.matrix(x)
49.  
    s <- as.matrix(subset)
50.  
    n <- nrow(x) ##n:样本量
51.  
    p <- ncol(x) ##p:变量数
52.  
     
53.  
    ##########变量选择矩阵的建立########
54.  
    if ((nrow(s)==1) && (ncol(s)==1)){ ##如果subset的形式是嵌套或不进行变量筛选（全模型平均）
55.  
    if (subset == 1){ ##如果是纯嵌套模型平均
56.  
    s <- matrix(1,nrow=p,ncol=p) ##生成（pxp）的矩阵，元素都是1 <-共p个备选模型
57.  
    s[upper.tri(s)] <- 0 ##矩阵的上三角部分改为0
58.  
    zero <- matrix(0,nrow=1,ncol=p) ##元素为0的（1xp）的矩阵
59.  
    s <- rbind(zero,s) ##将zero和s两个矩阵上下拼接 <-变量选择矩阵
60.  
    }
61.  
    if (subset == 2){ ##如果是全模型平均
62.  
    s <- matrix(0,nrow=2^p,ncol=p) ##生成一个（2^p x p）的矩阵，元素都是0 <-共2^p个备选模型
63.  
    s0 <- matrix(c(1,rep(0,p-1)),1,p) ##生成一个(1xp)的矩阵，[1,1]是1，后面全是0
64.  
    s1 <- matrix(c(rep(0,p)),1,p) ##生成一个(1xp)的矩阵，元素为0
65.  
    for (i in 2:2^p){ ##从2:2^p
66.  
    s1 <- s0 s1 ##初始化s1，在上一个s1的基础上加了一个s0，利用下面这个
67.  
    ####2进制推动变量选择矩阵的变化
68.  
    for (j in 1:p){ ##从1:p |
69.  
    if (s1[1,j] == 2){ ##如果s1[1,j]是2 |<- 2进制机制
70.  
    s1[1,j 1] <- s1[1,j 1] 1 ##s1[1,j 1]在原来的基础上加1 | 利用这个做变量选择
71.  
    s1[1,j] <- 0 ##s1[1,j]赋值0 |
72.  
    }
73.  
    }
74.  
    s[i,] <- s1 ##将s1存到s的第i行
75.  
    }
76.  
    }
77.  
    }
78.  
     
79.  
     
80.  
    #########
81.  
    m <- nrow(s) ##m是s的行数，备选模型的个数
82.  
    bbeta <- matrix(0,nrow=p,ncol=m) ##bbeta是（pxm）的矩阵，元素为0 <-所有备选模型的系数
83.  
    if (method == 2) ee <- matrix(0,nrow=n,ncol=m) ##如果是刀切模型平均，ee是（nxm）矩阵，元素为0
84.  
     
85.  
    for (j in 1:m){ ##从1:m,遍历每一个备选模型
86.  
    ss <- matrix(1,nrow=n,ncol=1) %*% s[j,] ##（nx1）元素为1的矩阵 X s[j,]
87.  
    indx1 <- which(ss[,]==1) ##indx1是被选中变量的索引
88.  
    xs <- as.matrix(x[indx1]) ##|
89.  
    xs <- matrix(xs,nrow=n,ncol=nrow(xs)/n) ##|xs是剔除未选中变量的样本矩阵
90.  
    if (sum(ss)==0){ ##如果备选模型中不含变量
91.  
    xs <- x ##样本矩阵xs
92.  
    betas <- matrix(0,nrow=p,ncol=1) ##系数矩阵为(px1)的0矩阵
93.  
    indx2 <- matrix(c(1:p),nrow=p,ncol=1) ##indx2是一个(px1)的矩阵，元素为1:p
94.  
    }
95.  
    if (sum(ss)>0){ ##如果备选模型含有至少一个变量
96.  
    betas <- solve(t(xs)%*%xs)%*%t(xs)%*%y ##第j个备选模型的估计系数
97.  
    indx2 <- as.matrix(which(s[j,]==1)) ##indx2是第j个候选模型包含的变量索引
98.  
    }
99.  
    beta0 <- matrix(0,nrow=p,ncol=1) ##beta0是(px1)的矩阵，元素为0
100.  
     beta0[indx2] <- betas ##|将第j个候选模型的系数填入bbeta的第j列
101.  
     bbeta[,j] <- beta0 ##|
102.  
     if (method == 2){ ##如果是刀切模型平均
103.  
     ei <- y - xs %*% betas ##ei是残差
104.  
     hi <- diag(xs %*% solve(t(xs) %*% xs) %*% t(xs))##hi是pm的n个对角元素
105.  
     ee[,j] <- ei*(1/(1-hi)) ##第j个备选模型的残差估计ee
106.  
     }
107.  
     }
108.  
      
109.  
     if (method == 1){ ##如果是MMA
110.  
     ee <- y %*% matrix(1,nrow=1,ncol=m) - x %*% bbeta ##ee是残差矩阵，(nxm)
111.  
     ehat <- y - x %*% bbeta[,m] ##第m个备选矩阵的残差估计
112.  
     sighat <- (t(ehat) %*% ehat)/(n-p) ##方差矩阵的估计值
113.  
     }
114.  
      
115.  
     a1 <- t(ee) %*% ee ##a1:(n-p)方差矩阵
116.  
     if (qr(a1)$rank<ncol(ee)) a1 <- a1 diag(m)*1e-10 ##由无穷小数将a1补成满秩的矩阵
117.  
     if (method == 1) a2 <- matrix(c(-sighat*rowSums(s)),m,1) ##MMA:a2:(mx1)矩阵，
118.  
     if (method == 2) a2 <- matrix(0,nrow=m,ncol=1) ##JMA:a2:(mx1)的0矩阵
119.  
     a3 <- t(rbind(matrix(1,nrow=1,ncol=m),diag(m),-diag(m))) ##
120.  
     a4 <- rbind(1,matrix(0,nrow=m,ncol=1),matrix(-1,nrow=m,ncol=1)) ##
121.  
      
122.  
     ###########最佳W的计算##########
123.  
     w0 <- matrix(1,nrow=m,ncol=1)/m ##w0是一个(mx1)的矩阵，元素为1/m
124.  
     QP <- solve.QP(a1,a2,a3,a4,1) ##计算argmin,利用quadprog包解决二次规划问题
125.  
     w <- QP$solution ##|将结果赋给w
126.  
     w <- as.matrix(w) ##|
127.  
     w <- w*(w>0) ##保留大于0的权重，小于0的记为0
128.  
     w <- w/sum(w0) ##?????sum(w0)不是1吗???
129.  
      
130.  
     ###########由最佳权重W计算得一堆结果########
131.  
     betahat <- bbeta %*% w ##估计系数
132.  
     ybar <- mean(y) ##y的均值
133.  
     yhat <- x %*% betahat ##y的估计值
134.  
     ehat <- y-yhat ##残差e的估计值
135.  
     r2 <- sum((yhat-ybar)^2)/sum((y-ybar)^2) ##R^2
136.  
     if (method == 1) cn=(t(w) %*% a1 %*% w - 2*t(a2) %*% w)/n ##MMA
137.  
     if (method == 2) cn=(t(w) %*% a1 %*% w)/n ##JMA
138.  
     list(betahat=betahat,w=w,yhat=yhat,ehat=ehat,r2=r2,cn=cn) ##显示结果
139.  
     }
140.  
      
141.  
      

 4.20更新

1.  
   ##HANSEN大佬的实验
2.  
   ##文章来源：LEAST SQUARES MODEL AVERAGING
3.  
   ##需要使用的gma函数先运行
4.  
    
5.  
   #用到的包
6.  
   library(MASS)
7.  
    
8.  
   #参数设定
9.  
   n <- 50 ##c(50,150,400,1000)
10.  
    alpha <- 0.5 ##c(0.5, 1.0, 1.0)
11.  
    M <- round(3*n^(1/3))
12.  
    p <- M
13.  
    beta <- matrix(1, nrow = p, ncol = 1)
14.  
    R2 <- seq(0.1,0.9,by=0.1)
15.  
    c <- ((1-R2)/R2)^(-0.5)
16.  
    t <- 1000000
17.  
    method=1;subset=1
18.  
    result <- matrix(nrow = t 1, ncol = 9)
19.  
    i=1
20.  
    a <- matrix(nrow = t 1, ncol = 9)
21.  
     
22.  
    options(digits = 8)
23.  
    for(T in 1:t){
24.  
    for (i in 1:9){
25.  
    for(j in 1:p){beta[j] <- c[i]*(2*alpha)^0.5*j^(-alpha-0.5)}
26.  
    x1 <- matrix(1, nrow = n, ncol = 1)
27.  
    x_sigma <- diag(1,p-1,p-1)
28.  
    x_mean <- matrix(0, nrow=p-1, ncol = 1)
29.  
    xp <- mvrnorm(n,x_mean,x_sigma)
30.  
    x <- cbind(x1,xp)
31.  
    e <- rnorm(n, 0 ,1)
32.  
    y <- x%*�ta e
33.  
    r <- gmaN(y,x,method,subset)
34.  
    d <- lm(y~x)$coefficients[2:11]
35.  
    d[is.na(d)] <- 0
36.  
    a[T,i] <- t(d-beta)%*%(d-beta)
37.  
    result[T,i]<-t(r$betahat-beta)%*%(r$betahat-beta)/a[T,i] ##<-这个地方有问题
38.  
    }
39.  
    }
40.  
     
41.  
    for (i in 1:9){
42.  
    result[t 1,i]<-mean(result[1:t,i])
43.  
    }
44.  
     
45.  
    plot(R2,result[t 1,],type = "l")
46.  
    points(R2,result[t 1,])

实验重复的是文章中的实验第一张图，但是文中说的评价标准是Risk(expected squared error)感觉是我理解有误，图片出来不对，后补。

但是流程是这么个流程没问题。

这篇好文章是转载于：学新通技术网