1.概率PCA建模

 概率PCA(Probability Principle Component Analysis)是从概率的角度理解主元分析这个事情。

 看到这里，我们会发现这个像是一个因子分析的模型，什么是因子分析，实际上就是一种数据简化技术，用较少的特征来反映原数据的基本结构。而这个较少的数据特征就是我们的隐变量，而原数据就是我们的。和相互独立，且都服从于高斯分布。而被称作因子载荷。这个模型意在用隐变量的特征来解释原数据的主要特征。
 但是这个和我们的降维有什么关系呢，我们如何通过这样一个模型来降维呢？换句话说我们的降维过程具体化到我们的目的是求什么？
 实际上，我们试图算一个后验,也就是得到一个条件分布，取到它的均值。要求后验，首先我们想到的是贝叶斯公式：

 由我们先前在(二)多元高斯分布与概率图条件独立性假设中提到的定理一可知仍然服从高斯分布：
因为：
;

因此：

的分布是我们的已知条件，那么剩下的是

因为：
所以:


我们令
 现在我们已经凑齐我们可以求解参数的要素了，我们需要显式的表达出这个后验概率，由于都服从于高斯分布，我们可以丢掉高斯分布前面的常数项:

 因为我们关心的已知的情况下的概率分布，所有与无关的项都视为常数，

 所以经过上面的凑项配方可以看出：

 因此:

 从这个高斯分布可以得出我们想要的最后结果:
 在这个式子里，哪些是我们知道的参数呢？数据是我们唯一知道的，都是我们需要估计的参数。接下来我们将用极大似然估计来讲这几个未知参数给估计出来。

2.极大似然估计求解

 我们想要用极大似然估计的话，我们就必须知道我们样本数据的概率分布，即 的分布，我们在上面求过这个分布：

 我们令
 似然函数为：

 为了求导方便，我们对似然函数取对数：

 我们先估计;

将这个结果带入对数似然函数中：

 由于后面一项是一个数，我们直接对数引入迹，对于迹来说

 这样操作之后，我们惊讶的发现

 我们得出的其中一项是协方差矩阵，所以可以得出我们的似然函数是：

 在对进行参数估计之前，我们可以先对C求导，因为C里同时包含了两个参数，通过对C和两个参数之间的关系，能直接得出两个参数的统计估计量：
 似然函数中涉及C的行列式，我们知道对行列式求导的公式为：

 所以对于原似然函数求导：

 我们的可以通过下式得来：
 因为：所以：

 代入原式中有：

 首先我们对进行参数估计：
 因为：

 所以：

 代入似然函数原式中：

 得出：

 即：

在进行下一步计算前，我们将给出下面等式并予以证明：

证明：
要证：

即证：

即证明：

上式恒成立，原式成立，证毕。
根据该等式，有：

 由于是的实对称矩阵，所以对其进行特征分解：

代入：

 这里我们可以看出WV还不是单位化的特征向量，因为：

 因此我们在两边都乘得：

 从上式中，我们得出是样本协方差矩阵的特征向量，而是它的特征值构成的对角矩阵。
 由于样本协方差矩阵也是一个实对称的半正定矩阵，我们也可以对其进行特征分解：

 解得

 所以：


 到这里我们得到的参数估计：

 但是这个式子里仍有未知参数,现估计:
 
 代入对求导的似然函数中：

 所以：
 根据Sherman-Morrison-Woodburg等式(这个博客开头有个地方和符号反了)

 利用上述不等式，我们继续计算：

 由于和都是实对称矩阵:

 由于上面我们已经得到：所以：
 将这个结论代入

 所以:

 因为:
 因为:
 即：

 至此，也估计出来了，那么估计的:

有了这些：

剩下就是代入的问题了。
展示代码如下：

3.代码

# #            ** Statistical Learning **
# #            File Name:     Probability Principle Component Analysis
# #            Author:        Pan
# #            Description:   a kind of dimension reduction based on Factor Analysis.
# #            Date:          2020/8/3
# #            **      --Start--      **
import numpy as np
class PPCA:
    def __init__(self,X,q):
        self.X=X
        self.mu=np.array([np.average(self.X,axis=0)]).T
        self.n=np.shape(self.X)[0]
        self.p=np.shape(self.X)[1]
        self.q=q
        self.H=np.identity(self.n)-(1/self.n)*np.ones([self.n,1]).dot(np.ones([self.n,1]).T)
        self.S=self.X.T.dot(self.H).dot(self.X)
    def dimensionality_reduction(self):
        self.eigenvalues,self.eigenvectors=np.linalg.eig(self.S)
        self.ranking=np.argsort(self.eigenvalues)
        self.eigenvalues_q=self.eigenvalues[np.where(self.ranking<self.q)]
        self.eigenvectors_q=self.eigenvectors[np.where(self.ranking<self.q)]
        self.Gamma=np.diag(self.eigenvalues)
        self.Gamma_q=np.diag(self.eigenvalues_q)
        self.sigma_square=(1/(self.p-self.q))*(np.sum(self.eigenvectors)-np.sum(self.eigenvectors_q))
        self.Lambda=self.Gamma_q-self.sigma_square*np.identity(self.q)
        self.W=self.eigenvectors_q.T.dot(self.Lambda**(1/2))
        self.Z=np.linalg.inv(self.W.T.dot(self.W)+self.sigma_square*np.identity(self.q)).dot(self.W.T).dot(self.X.T-self.mu)
        return self.Z.T
if __name__ == '__main__':
    import seaborn as sns
    import matplotlib.pyplot as plt
    from sklearn.datasets import load_diabetes
    def plot(dfData):
        plt.subplots(figsize=(9, 9))
        sns.heatmap(dfData, annot=True, vmax=1, square=True, cmap="Blues")
        plt.show()
    data=load_diabetes()
    X=data["data"]
    PPCA_Algorithm=PPCA(X,4)
    Z=PPCA_Algorithm.dimensionality_reduction()
    plot(np.corrcoef(X.T))
    plot(np.corrcoef(Z.T))

结果展示
1.降维前：

降维前
2.降维后：
降维后

4.核概率PCA

 因为，这里忘记说了，之所以可以取单位矩阵是因为，V是特征分解出来的，而W是未知的，只要求V是正交矩阵都行，为了方便，直接就取单位矩阵。这样我们就可以得出:
 
 令:
 将进行特征分解可得

 从中可以得出就是协方差矩阵特征向量构成的正交矩阵，即：

 引入核:其中映射g属于再生核希尔伯特空间。

 将结果带入概率PCA的最后公式：

上代码：

import numpy as np
class PPCA:
    def __init__(self,X,q):
        self.X=X
        self.mu=np.array([np.average(self.X,axis=0)]).T
        self.n=np.shape(self.X)[0]
        self.p=np.shape(self.X)[1]
        self.q=q
        self.H=np.identity(self.n)-(1/self.n)*np.ones([self.n,1]).dot(np.ones([self.n,1]).T)
        self.S=self.X.T.dot(self.H).dot(self.X)
    def Guass_Kernal(self,sigma):
        self.sigma=sigma
        E_minus_w_square = (np.repeat(np.exp(-1 * (np.sum(self.W.T * self.W.T, axis=1))), self.q) / self.sigma).reshape(
            self.q, self.q)
        E_w_wT = np.exp(self.W.T.dot(self.W) / self.sigma)
        return E_minus_w_square*E_w_wT*E_minus_w_square.T
    def Multinomial_Kernel(self,M_n,alpha):
        return (self.W.T.dot(self.W)+alpha)**M_n
    def dimensionality_reduction(self,):
        self.eigenvalues,self.eigenvectors=np.linalg.eig(self.S)
        self.ranking=np.argsort(self.eigenvalues)
        self.eigenvalues_q=self.eigenvalues[np.where(self.ranking<self.q)]
        self.eigenvectors_q=self.eigenvectors[np.where(self.ranking<self.q)]
        self.Gamma=np.diag(self.eigenvalues)
        self.Gamma_q=np.diag(self.eigenvalues_q)
        self.sigma_square=(1/(self.p-self.q))*(np.sum(self.eigenvectors)-np.sum(self.eigenvectors_q))
        self.Lambda=self.Gamma_q-self.sigma_square*np.identity(self.q)
        self.inv_lambda=np.linalg.inv(self.Lambda)
        self.W=self.eigenvectors_q.T.dot(self.Lambda**(1/2))
        self.T=(self.Gamma_q**(1/2)).dot(self.inv_lambda**(1/2)).dot(self.Guass_Kernal(0.01)).dot(self.inv_lambda**(1/2)).dot(self.Gamma_q**(1/2))
        _, self.M=np.linalg.eig(self.T)
        self.Z=self.H.dot(self.X).dot(self.eigenvectors_q.T).dot(self.M.T).dot(np.linalg.inv(self.Gamma_q)**(1/2)).dot(self.Lambda**(1/2)).dot(np.linalg.inv(self.Guass_Kernal(0.01)+self.sigma_square*np.identity(self.q)))
        return self.Z
if __name__ == '__main__':
    import seaborn as sns
    import matplotlib.pyplot as plt
    from sklearn.datasets import load_diabetes
    def plot(dfData):
        plt.subplots(figsize=(9, 9))
        sns.heatmap(dfData, annot=True, vmax=1, square=True, cmap="Reds")
        plt.show()
    data=load_diabetes()
    X=data["data"]
    print(X.shape)
    PPCA_Algorithm=PPCA(X,4)
    Z=PPCA_Algorithm.dimensionality_reduction()
    plot(np.corrcoef(X.T))
    plot(np.corrcoef(Z.T))

事实上，我使用了高斯核和多项式核，发现与没有引入之前相比，在数据集diabetes上并没有多大差别。有时候特征分解还有复数出现，不知道如何处理，希望以后能解决吧。
结果：

KPPCA.png

（OVER）