切片逆回归

简介

切片逆回归(Slice Inverse Regression, SIR)是Li(1991)¹提出的一种经典的充分降维方法。SIR探究了多维变量 $x$ 对单变量 $y$ 的回归而不是一维变量 $y$ 对多维变量 $x$ 的回归。在对多维解释变量 $x$ 的降维过程中充分利用了 $y$ 的信息。

一般模型

一般模型： $y=f(\beta^{T}_{1}x,\dots,\beta_{K}^{T}x,\varepsilon)$
- $\varepsilon$ 是一个独立于 $x$ 的随机变量， $f$ 为一个未知的连接函数；
- $y$ 关于 $x$ 的条件分布可以由 $x$ 的 $K$ 个线性组合 $\beta^{T}_{1}x,\dots,\beta_{K}^{T}x$ 得到，并不会损失 $x$ 中包含的原始信息；
- 等价于在给定 $x$ 的 $K$ 个线性组合时，响应变量 $y$ 与解释变量 $x$ 是独立的。当 $K$ 远小于 $x$ 的维度 $p$ 时，便达到降维目的;
基本定理：
- 逆回归曲线 $E (X ∣ Y) ? E (X)$ 被包含在 $\Sigma_{x}\beta_{i},i=1\,\dots,k$ 张成的线性空间中；
- 在线性条件下，满足 $\Sigma_{\eta}b_{i}=\lambda_{i}\Sigma_{x}b_{i}$ 的非零相对特征 $\lambda_{i}$ 不超过 $k$ 个，其对应的相对特征向量 $b_{i}$ 被包含在 $\beta_{i},i=1,\dots,k$ 张成的线性空间中；

具体步骤

切片逆回归对数据 $(Y_{i},X_{i}),(i=1,\dots,n)$ 进行操作;
通过仿射变换标准化 $X_{i}$ 得到 $\widetilde{X_{i}}=\widehat{\Sigma}_{xx}^{-1/2}(X_{i}-\overline{X}),i=1,\dots,n$ ，其中 $\widehat{\Sigma}_{xx}$ 和 $\overline{X}$ 分别是 $X$ 的样本协方差矩阵和样本均值;
将 $Y$ 的取值 $Y_{i}$ 从小到大进行排序，并切为 $H$ 片，记为 $I_{1},I_{2},\dots,I_{h}$ 。其中 $I_{h}$ 包含 $Y_{i}$ 的概率为 $P_{h}$ ：
$P_{h}=\frac{1}{n}\sum_{i=1}^{n}\delta_{h}(Y_{i})$

$\delta_{h}\left(Y_{i}\right)=\left\{\begin{array}{l} 0, \text { if } Y_{i} \text { inside the } \mathrm{h} \\ 1, \text { if } Y_{i} \text { outside the } \mathrm{h} \end{array}\right.$
对于每一个切片，计算其 $x_{i}$ 的样本均值 $m_{h}$ ，即 $m_{h}=\frac{1}{nP_{h}}\sum_{y\in{I{_h}}}\widetilde{X_{i}}$ ;
对数据 $m_{h}$ 进行（加权）主成分分析，加权的协方差矩阵为 $V=\sum_{h=1}^{H}p_{h}m_{h}m_{h}^{'}$ ，并找出V的特征值和特征向量。
令 $k$ 个最大的特征向量为 $\eta_{k},(k=1,2,\dots,k)$ ，可得 $\beta_{k}=\eta_{k}\Sigma_{xx}^{-\frac{1}{2}},k=1,2,\dots,k$

切片逆回归的Python实现

仅包含一般情况，稀疏矩阵代补充

class SIR:
    def __init__(self, H, K, estdim=0, Cn=1):
        self.H = H
        self.K = K
        self.estdim = estdim
        self.Cn = Cn

    def fit(self, X, Y):
        self.X = X
        self.Y = Y
        self.mean_x = np.mean(X, axis=0)  # 计算样本均值
        self.sigma_x = np.cov(X, rowvar=False)  # 计算样本协方差矩阵
        self.Z = np.matmul(X - np.tile(self.mean_x, (X.shape[0], 1)),
                           fractional_matrix_power(self.sigma_x, -0.5))  # 标准化后的数据阵
        n, p = self.Z.shape
        if self.Y.ndim == 1:  # 判断响应变量Y的维度
            self.Y = self.Y.reshape(-1, 1)
        ny, py = self.Y.shape
        W = np.ones((n, 1)) / n
        nw, pw = W.shape
        # 输入数据异常判断
        if n != ny:
            raise ValueError('X and Y must have the same number of samples')
        elif p == 1:
            raise ValueError('X must have at least 2 dimensions')
        elif py != 1:
            raise ValueError('Y must have only 1 dimension')
        # 将y分为H片，c为每个片的样本数
        c = np.ones((1, self.H)) * (n // self.H) + np.hstack(
            [np.ones(
                (1, n % self.H)), np.zeros((1, self.H - n % self.H))])
        cumc = np.cumsum(c)  # 计算切片累计和
        # 参照Y的取值从小到大进行排序
        temp = np.hstack((self.Z, self.Y, W))
        temp = temp[np.argsort(temp[:, p])]
        # 提取排序后的z,y,w
        z, y, w = temp[:, :p], temp[:, p:p + 1], temp[:, p + 1]
        muh = np.zeros((self.H, p))
        wh = np.zeros((self.H, 1))  # 每个切片的权重
        k = 0  # 初始化切片编号
        for i in range(n):
            if i >= cumc[k]:  # 如果超过了切片的边界，则换下一个切片
                k += 1
            muh[k, :] = muh[k, :] + z[i, :]  # 计算切片内自变量之和
            wh[k] = wh[k] + w[i]  # 计算每个切片包含Yi的概率
        # 计算每个切片的样本均值,将其作为切片内自变量的取值
        muh = muh / (np.tile(wh, (1, p)) * n)
        # 加权主成分分析
        self.M = np.zeros((p, p))  # 初始化切片 加权协方差矩阵
        for i in range(self.H):
            self.M = self.M + wh[i] * muh[i, :].reshape(-1, 1) * muh[i, :]
        if self.estdim == 0: # 一般情况
            self.D, self.V = eigs(A=self.M, k=self.K, which='LM')
        else:
            """ # 稀疏矩阵情况，待修正
            [V D] = np.linalg,eig(full(M))
            lambda = np.sort(abs(diag(D)),'descend')
            L = np.log(lambda+1) - lambda
            G = np.zeros((p,1))
            if Cn == 1
                Cn = n^(1/2)
            elif Cn == 2
            Cn = n^(1/3)
            elif Cn == 3:
                Cn = 0.5 * n^0.25
            for k in range(p):
                G(k) = n / 2 * sum(L(1:k)) / sum(L) - Cn * k * (k+1) / p
            maxG, K = np. max(G)
            V, D = eigs(M,K,'lm')
            """
            pass
        return self.V, self.K, self.M, self.D
    def transform(self):
        hatbeta = np.matmul(fractional_matrix_power(self.sigma_x, -0.5),self.V)
        return hatbeta

X = pd.read_csv('X.csv',header=None).values
Y = pd.read_csv('Y.csv',header=None).values
model = SIR(H=6, K=4, estdim=0, Cn=1)
V, K, M, D = model.fit(X,Y)
hatbeta = model.transform()

Li, Ker-Chau. “Sliced Inverse Regression for Dimension Reduction: Rejoinder.” Journal of the American Statistical Association, vol. 86, no. 414, 1991, pp. 337–42. JSTOR, https://doi.org/10.2307/2290568. Accessed 17 Jul. 2022. ??