精华内容
参与话题
问答
  • 函数型数据主成分分析

    千次阅读 2019-04-15 08:59:29
    函数型主成分分析(FPCA,Functional Principal Components Analysis)是传统的PCA的一种推广。 考虑我们已经从数据中得到拟合曲线xi(s),s∈T,i=1,⋯ ,nx_{i}(s), s \in \mathcal{T}, i=1, \cdots, ...

    函数型数据主成分分析

    基本思想

    函数型主成分分析(FPCA,Functional Principal Components Analysis)是传统的PCA的一种推广。

    考虑我们已经从数据中得到拟合曲线xi(s),sT,i=1, ,nx_{i}(s), s \in \mathcal{T}, i=1, \cdots, n,所谓的第一主成分,就是我们希望能找到一个模为1的函数β(s)\beta(s),使得{xi}\{x_i\}β\beta上的投影(L2L_2內积){ξi}\{\xi _i\}的方差达到最大,方差最大其实也就体现{xi}\{x_i\}整体到β\beta的距离达到最小。β\beta一般就叫做权重函数(可以理解为“坐标轴”单位长度量)。

    我们管各个函数到β\beta上的投影叫做观测曲线的主成分得分:
    ξi=Tβ(s)xi(s)ds,i=1, ,n \xi_{i}=\int_{\mathcal{T}} \beta(s) x_{i}(s) d s, \quad i=1, \cdots, n

    故而,求解第一个主成分就变成了求解一个优化问题:
    max1ni=1nξi2=max1ni=1n(Tβ(s)xi(s)ds)2 s.t. β2=Tβ(s)β(s)ds=1 \begin{aligned} \max \frac{1}{n} \sum_{i=1}^{n} \xi_{i}^{2} &=\max \frac{1}{n} \sum_{i=1}^{n}\left(\int_{\mathcal{T}} \beta(s) x_{i}(s) d s\right)^{2} \\ \text { s.t. } &\|\beta\|^{2}=\int_{T} \beta(s) \beta(s) d s=1 \end{aligned}
    求解这个优化问题,我们就得到了第一主成分β1(s)\beta^1(s)

    kk主成分无非就是在满足和前面k1k-1个主成分权重函数垂直的基础上,求解上述优化问题而已,即求解

    max1ni=1nξi2=max1ni=1n(Tβ(s)xi(s)ds)2 s.t. β2=Tβ(s)β(s)ds=1Tβ(s)βl(s)ds=0,l=1, ,k1 \begin{array}{l}{\max \frac{1}{n} \sum_{i=1}^{n} \xi_{i}^{2}=\max \frac{1}{n} \sum_{i=1}^{n}\left(\int_{\mathcal{T}} \beta(s) x_{i}(s) d s\right)^{2}} \\ {\text { s.t. }\|\beta\|^{2}=\int_{T} \beta(s) \beta(s) d s=1} \\ {\int_{T} \beta(s) \beta^{l}(s) d s=0, l=1, \cdots, k-1}\end{array}

    这个优化问题的解可以表述如下。记协方差函数:

    v(s,t)=1n1i=1n(xi(s)x(s))(xi(t)x(t)) v(s, t)=\frac{1}{n-1} \sum_{i=1}^{n}\left(x_{i}(s)-\overline{x}(s)\right)\left(x_{i}(t)-\overline{x}(t)\right)

    那么权重函数满足特征方程:

    Tv(s,t)β(t)dt=λβ(s) \int_{\mathcal{T}} v(s, t) \beta(t) d t=\lambda \beta(s)

    定义积分变换:

    Vβ(s)=Tv(s,t)β(t)dt V \beta(s)=\int_{\mathcal{T}} v(s, t) \beta(t) d t

    这里的VV称为协方差算子,它将函数β\beta变成一个函数。那么,我们有:

    Vβ(s)=λβ(s) V \beta(s)=\lambda \beta(s)

    我们也类比PCA,使用特征值的累积贡献率来衡量主成分所占比例:

    FVE=i=1Kλi/i=1n1λi \mathrm{FVE}=\sum_{i=1}^{K} \lambda_{i} / \sum_{i=1}^{n-1} \lambda_{i}

    这里之所以对λ\lambda只累计到nn是因为协方差算子VV的秩为样本数量减一个,则非零特征根的个数最多为n1n-1个。

    特征问题求解

    由上述已知,我们求解主成分最后归结为求解一个特征值问题。

    求解这个问题,目前比较流行的有三种方法:

    • 对函数进行SVD离散化
    • 对函数进行基函数展开
    • 运用一般性的数值积分方法

    我们最后需要的是特征函数,为了避免插值而带来更大的误差,我选用对基函数进行展开的方法。下面简单介绍一个对函数进行基函数展开的基本思路。

    我们的样本基函数xix_i可以通过基函数展开,如下:
    Xi(s)=k=1KcikΦk(s),i=1,2,,N X_{i}(s)=\sum_{k=1}^{K} c_{i k} \Phi_{k}(s), i=1,2, \ldots, N

    我们记X=(x1,x2,,xN),Φ=(Φ1,,Φk),C=(cik)N×K X=\left(x_{1}, x_{2}, \ldots, x_{N}\right)^{\prime}, \Phi=\left(\Phi_{1}, \ldots, \Phi_{k}\right)^{\prime}, C=\left(c_{i k}\right)_{N \times K}

    那么样本函数就可以写为等价的矩阵形式X=CΦX=C \Phi。那么协方差函数就可以写为(假设已经标准化):

    v(s,t)=1n1Φ(s)CCΦ(t) v(s, t)=\frac{1}{n-1} \Phi^{\prime}(s) C^{\prime} C \Phi(t)

    定义K阶对称矩阵

    W=ΦΦ W=\int \Phi \Phi^{\prime}

    当选择正交基的时候,比如说正交傅里叶基,这就是一个单位矩阵。关于这个基如何选取,我们后面还会详谈。

    同样地,将特征函数进行展开:

    β(s)=k=1KbkΦk(s)=Φ(s)b \beta(s)=\sum_{k=1}^{K} b_{k} \Phi_{k}(s)=\Phi^{\prime}(s) b

    将其代入

    Tv(s,t)β(t)dt=λβ(s) \int_{\mathcal{T}} v(s, t) \beta(t) d t =\lambda \beta(s)

    就可以得到(N=n1N=n-1):

    1NΦ(s)CCWb=λΦ(s)b \frac{1}{N} \Phi^{\prime}(s) C^{\prime} C W b=\lambda \Phi^{\prime}(s) b

    进一步能得到1NCCWb=λb\frac{1}{N} C^{\prime} C W b=\lambda b,由特征向量正交和单位长度的约束要求,有bkWbk=1,bkWbm=0,kmb_{k}^{\prime} W b_{k}=1, b_{k}^{\prime} W b_{m}=0,k \neq m

    WW做cholesky分解,可得W=LLW=LL'

    定义u=Lbu=L'b,那么上述问题就变成了对称矩阵的代数特征值问题:

    1NLCCLu=λu \frac{1}{N} L' C^{\prime} C L u=\lambda u

    据此可以求得uu,进而求得bb,最后求得特征函数β\beta

    关于基函数的选取

    如上所述,可以看出,如果所选的基函数是正交的,本质上和PCA的以拟合系数为坐标点的函数空间PCA推广是实际上是一样的。若基函数不是正交的,无非就是在此基础上对要求特征值的矩阵得多乘一个W=ΦΦW=\int \Phi \Phi^{\prime},再求特征向量,以及进行WW意义下对特征向量进行单位化而已(不单位话也没事,只不过权重函数β\beta不再是模长为1的而已)。

    常用的基函数有傅里叶基函数和B样条基函数,傅里叶基函数适用于周期性函数数据,B样条基函数适用于非周期函数数据,当然,也可以用多项式基函数。

    B样条基函数的递归定义为:

    Bj,0(x)={1,tjx<tj+10,elseBi,k(x)=xtiti+ktiBi,k1(x)+ti+k+1xti+k+1ti+1Bi+1,k1(x),k>0 \begin{array}{c}{B_{j, 0}(x)=\left\{\begin{array}{l}{1, t_{j} \leq x<t_{j+1}} \\ {0, \text {else}}\end{array}\right.} \\ {B_{i, k}(x)=\frac{x-t_{i}}{t_{i+k}-t_{i}} B_{i, k-1}(x)+\frac{t_{i+k+1}-x}{t_{i+k+1}-t_{i+1}} B_{i+1, k-1}(x), k>0}\end{array}

    程序代码

    下面附一段简单的以多项式为基的matlab代码。

    clc
    clear
    close all
    format short
    digits(5)
    %% 设定维数
    p = 8;
    d = 2;
    %% 读入数据,保存
    D = dir('../data/*.csv');
    D = struct2cell(D);%结构体不好调用,转成元胞数组来使用
    csvs_name = D(1,:);
    N = 0;%这个表示数据的个数
    %csvs_name = csvs_name(5:6);
    for csv_name = csvs_name %cpicName = 'a001.bmp'
        N = N+1;
        csvName = csv_name{1};
        Data{N} = load(['../data/' csvName]);
    end
    %clear csv_name csvN csvName csvs_name D
    %% 定义原空间基,生成模型函数,由模型函数通过拟合得到点的坐标表示
    %p = 15;%原空间维数
    %%根据p自动生成模型函数
    eval_poly = 'a(1)';
    for i = 2:p
        eval_poly = [eval_poly '+a(' num2str(i) ').*x.^' num2str(i-1)];
    end
    eval_poly = ['modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(' eval_poly ');'];
    eval(eval_poly);
    
    %modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*sin(x)+a(3).*cos(x)+a(4).*1./x+a(5).*x+a(6).*x.^2);
    %   modelfun = @(a,x)(x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*x+a(3).*x.^2+...
    %       a(4).*x.^3))%+a(5).*x.^4+a(6)*x.^5))+...
    %       a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
    %         modelfun = @(a,x)((a(1)+a(2).*x.^2+a(3).*x.^2+...
    %         a(4).*x.^3+a(5).*x.^4+a(6)*x.^5+...
    %         a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
    
    %%做一波拟合,求出系数,进行拟合,拟合结果为A
    %%A中存的是数据点的坐标表示
    for k = 1:N
        data = Data{k};
        data(:,2) = abs(data(:,2));%取绝对值
        x = data(:,1);
        y = data(:,2);
        %y = y./(x.*(1-x));%%%%%%%%如果将x(1-x)项挪到左边,拟合会出问题
        a0 = ones(1,p);
        a = nlinfit(x,y,modelfun,a0);
        A(:,k) = a';
    end
    %clearvars -except A modelfun x y
    
    %% 开始奇异值分解降维,得到新空间的原点和坐标轴,坐标表示
    %%由奇异值分解做主成分分析,A的每一列为原空间下的系数,U的每一列为主成分
    [p,N] = size(A);
    %d = 5;%新空间维数
    x0 = mean(A,2);%计算样本均值
    A_shift = A - repmat(x0,1,N);%中心平移每个训练样本,A_shift为原数据点的平移表示
    
    %%%%对于A_shift的补充处理
    syms x;
    base_vec_T = x.^(0.5).*(1-x).^(0.5).*x.^[0:p-1];
    base_vec = conj((base_vec_T))';
    W = int(base_vec*base_vec_T,[0,1]);
    L_prime = chol(W);
    L_prime_A_shift = L_prime*A_shift;
    
    %%%%
    [U,S,V] = svd(L_prime_A_shift);
    B = L_prime\U;
    if size(S,2) == 1
        S_diag = S(1);
    else
        S_diag = diag(S);
    end
    R = cumsum(S_diag)./sum(S_diag);
    Coord_new = B(:,1:d);
    
    
    %% 新空间的坐标轴和原点的函数表示
    %%生成新的空间中的基函数,找到d个新空间基函数及原点坐标
    %%Coord_new x0中存的是新坐标系,fi存的是新的函数空间的原点和坐标轴
    for k = 1:d;%取出第k个基函数
        syms x;
        c = B(:,k);
        exp = modelfun(c,x);
        exp = simplify(exp);%第k个基函数的表达式
        f{k} = exp;
        digits(2);
        latex(vpa(exp,2));%写成latex表达式黏贴到解释器中
    end
    
    f0 = modelfun(x0,x);
    f0 = simplify(f0);
    latex(vpa(f0,2));
    % for i = 1:d
    %     f_fun{i} = matlabFunction(f{i});
    % end
    % f0_fun = matlabFunction(f0);
    
    %% 求原坐标系中的投影后在原空间中的位置
    
    
    
    %A_new_coord = Coord_new'*A_shift;%A_new_coord表示原数据点在新空间中的表示
    
    syms x;
    for k = 1:N
        a = A_shift(:,k);
        pk_f = modelfun(a,x);
        for j = 1:d
            A_new_coord(j,k) = int(pk_f.*f{j},0,1);
        end
    end
    
    %%函数运算
    for k = 1:N
        A_shadow_f = f0;
        for i = 1:d
            A_shadow_f = A_shadow_f+f{i}*A_new_coord(i,k);
        end
        A_shadow_f = simplify(A_shadow_f);
        A_shadow_fs{k} = A_shadow_f;
        A_shadow_fs_fun{k} = matlabFunction(A_shadow_f);
    end
    
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % %%坐标运算,可以注释掉
    % A_shadow_coord = Coord_new*A_new_coord+repmat(x0,1,N);%A_shadow_coord表示投影坐标点在原坐标系下的表示
    % syms x;
    % for k = 1:N
    %     a = A_shadow_coord(:,k);
    %     A_shadow_coord_f = modelfun(a,x);
    %     A_shadow_coord_fs{k} = simplify(A_shadow_coord_f);
    % end
    %
    % %%不考虑舍入误差的情况下,可以证明坐标运算和函数运算的得到的函数表达式是一样的,验证代码如下
    % k=5;
    % diff = A_shadow_coord_fs{k}-A_shadow_fs{k}
    % diff = simplify(diff)
    % digits(2);
    % latex(vpa(diff,2))
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %% 画图和求误差
    Errors = [];
    NN = 1000;
    h = 1/NN;
    xx = linspace(0,1,NN+1);
    close all;
    for k = 1:N
        %%原数据点
        data = Data{k};
        data(:,2) = abs(data(:,2));%取绝对值
        x = data(:,1);
        y = data(:,2);
        %%拟合曲线
        a = A(:,k);
        data_fit_f = @(x) modelfun(a,x);
        %%投影函数
        A_shadow_f_fun = A_shadow_fs_fun{k};
        %%原点函数
        f0_fun = matlabFunction(f0);
        
        %%开始画图
        figure(k);
        scatter(x,y,'.','MarkerEdgeColor','b',...
                  'MarkerFaceColor','c',...
                  'LineWidth',0.5);
        hold on;
        
        y_fit = data_fit_f(xx);
        y_shadow = A_shadow_f_fun(xx);
        y_f0 = f0_fun(xx);
        plot(xx,y_fit,'g');
        plot(xx,y_shadow,'black');
        plot(xx,y_f0,'r');
        legend('数据','拟合','投影','原点');
        
        %%求拟合曲线和投影函数之间的L2误差
        %Errors(end+1) = norm((y_fit-y_shadow),2);
        Errors(end+1) = sum(abs(y_fit-y_shadow)*h);
    end
    
    Errors_mean = mean(Errors)
    
    
    
    %% test
    % close all;
    % plot(xx,xx);
    % hold on;
    % plot(xx,xx.^9)
    
    展开全文
  • 函数型数据主成分分析(FPCA)

    千次阅读 2019-11-24 19:57:27
    更近一步,我们介绍了函数型主成分分析方法(FPCA),包括其基本思想、数学推导、算法描述等,最为重要的是,我们将该方法和本领域进行结合,有了一些新的思考,感谢"数据科学与矩阵优化"课程给带来的灵感。...

    本文主要介绍了以下几个方面的内容:简单介绍了经典的主成分分析方法,包括其数学推导,算法步骤,和几个实际算例;简单介绍了其它的数据降维方法,譬如局部线性嵌入以及它的简单算例;更近一步,我们介绍了函数型主成分分析方法(FPCA),包括其基本思想、数学推导、算法描述等,最为重要的是,我们将该方法和本领域进行结合,有了一些新的思考。

    前言

    “维数灾难"带来的直接结果就是很多低维空间行之有效的算法在高维空间中变得不可计算,为此,我们需要进行降维。在另一个方面,数据偏平化的情况下,降维有助于我们抓住数据的主要结构,过滤可能的误差带来的影响,使模型更加真实。另外,在某些情况下,降维可用于可视化。数据降维的方法有很多,比如说基于"最小化投影误差”(最大化类内方法)的主成分分析方法(PCA),以及基于保持拓扑结构不变(高维空间中是邻居,到了地位空间中还是邻居)的局部线性嵌入(LLE)等方法。

    在多元统计分析中,主成分分析(Principal Components Analysis,PCA)是一种统计分析、简化数据集的方法。它利用正交变换来对一系列可能相关的变量的观测值进行线性变换,从而投影为一系列线性不相关变量的值,这些不相关变量称为主成分(Principal Components)。具体地,主成分可以看做一个线性方程,其包含一系列线性系数来指示投影方向。PCA对原始数据的正则化或预处理敏感(相对缩放)。PCA本质上寻求的是数据点在低秩空间中的一个表示。

    高维数据,意味着数据需要多于两个或三个维度来表示,一般很难被解释。一种简化的方法是假设数据嵌在高维空间的一个非线性流形上。如果这个流形维数足够低,那么数据可以在低维空间中被可视化。局部线性嵌入(Locally Linear Embedding,LLE)关注于降维时保持样本局部的线性特征,由于LLE在降维时保持了样本的局部特征,它广泛的用于图像图像识别,高维数据可视化等领域。

    有了降维和主成分分析,我们做PDE的就会思考,既然可以对Rn\mathbb{R}^n空间中的数据做降维,那么函数作为一组基函数的线性组合,如果将基函数看作一些坐标系中的一个个坐标轴,是否也可以对函数空间中的"数据"做降维呢?答案是肯定的。函数型(数据)主成分分析(Functional Principal Components analysis)可以视为是传统的主成分分析的一种推广。类比于PCA,它希望能将高维函数空间中的函数放到低维空间中去表示,而使得被表示的数据集损失最小。更通俗地说,就是希望用更少的基函数来表示某个已知基函数的函数空间的一堆函数,新空间的基函数用旧空间的基函数来线性表出。那么,我们就需要定义函数之间的距离,函数空间的內积等等。

    主成分分析(PCA)

    数据降维简介

    在机器学习和统计学领域,降维是指在某些限定条件下,降低随机变量个数,得到一组"不相关"主变量的过程。
    降维可进一步细分为变量选择和特征提取两大方法。除了考虑"维数灾难"的问题,降维还有一些本质的原因。目前大部分降维算法处理向量表达的数据,也有一些降维算法处理高阶张量表达的数据。之所以使用降维后的数据表示是因为在原始的高维空间中,包含有冗余信息以及噪音信息,在实际应用例如图像识别中造成了误差,降低了准确率,通过降维,我们希望减少冗余信息所造成的误差,提高识别(或其他应用)的精度,又或者希望通过降维算法来寻找数据内部的本质结构特征。目前比较流行的降维算法有主成分分析、线性判别分析、局部线性嵌入和拉普拉斯特征映射等等。

    PCA算法的原理解释

    所谓的主成分分析,不过是在高维的空间中寻找一个低维的正交坐标系,比如说在三维空间中寻找一个二维的直角坐标系。那么这个二维的直角坐标系就会构成一个平面,将三维空间中的各个点在这个二维平面上做投影,就得到了各个点在二维空间中的一个表示,由此数据点就从三维降成了二维。

    这个过程的关键在于,我们如何选取这个低维的坐标系,即坐标系原点在哪?各个轴朝哪个方向?一个原则就是使得各个点到到这个平面的距离平方和达到最小。由此,通过简单地数学推导,就能得到原点的表达公式和坐标系的各个基向量的表达公式。

    PCA算法的数学推导

    我们假设输入为p维的N个对象,XX表示如下图所示的一个矩阵:
    在这里插入图片描述
    通过PCA降维,将其降为d维的N个对象,假设为YY,同前,每列表示一个对象,每行表示一个特征:
    在这里插入图片描述

    我们要将所有点投影到新的坐标系中去,无非是寻找新坐标系的坐标原点和各个坐标轴。

    我们假设WW的每一列为新的坐标系中单位正交的坐标轴表示,x0x_0为新坐标系的原点(相对于原坐标系)。
    那么,我们要做的就是找到一个合适的WWx0x_0,使其极小化所有点到新的坐标平面的距离平方和。
    容易知道,每一个点到新坐标系的距离平方为(其中X=(x0,W)\underline X = (x_0,W)表示的是位置参数):
    DsX(x,X)=(xx0i=1dwiT(xx0)wi)T(xx0i=1dwiT(xx0)wi)\operatorname{Ds}_{X}(x, \underline{X})=\left(x-x_{0}-\sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right) w_{i}\right)^{T}\left(x-x_{0}-\sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right) w_{i}\right)
    对其进行化简,可得:
    DsX(x,X)=(xx0i=1dwiT(xx0)wi)T(xx0i=1dwiT(xx0)wi)=(xx0)T(xx0)(i=1dwiT(xx0)wi)T(xx0)(xx0)Ti=1dwiT(xx0)wi+(i=1dwiT(xx0)wi)Ti=1dwiT(xx0)wi\begin{array}{l}{D s_{X}(x, \underline{X})=\left(x-x_{0}-\sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right) w_{i}\right)^{T}\left(x-x_{0}-\sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right) w_{i}\right)} \\ {=\left(x-x_{0}\right)^{T}\left(x-x_{0}\right)-\left(\sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right) w_{i}\right)^{T}\left(x-x_{0}\right)} \\ {-\left(x-x_{0}\right)^{T} \sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right) w_{i}+\left(\sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right) w_{i}\right)^{T} \sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right) w_{i}}\end{array}
    进而有,
    DsX(x,X)=(xx0)T(xx0)i=1dwiT(xx0)(xx0)Twi\operatorname{Ds}_{X}(x, \underline{X})=\left(x-x_{0}\right)^{T}\left(x-x_{0}\right)-\sum_{i=1}^{d} w_{i}^{T}\left(x-x_{0}\right)\left(x-x_{0}\right)^{T} w_{i}
    让所有点到投影点距离平方和最小,即求解约束优化问题:
    minXkDsX(xk,X)=k(xkx0)T(xkx0)i=1dwiTk(xkx0)(xkx0)TwiwiTwj=δijδij=1,i=jδij=0,ij\begin{array}{l}{\min _{\underline{X}} \sum_{k} D s_{X}\left(x_{k}, \underline{X}\right)=\sum_{k}\left(x_{k}-x_{0}\right)^{T}\left(x_{k}-x_{0}\right)} \\ {-\sum_{i=1}^{d} w_{i}^{T} \sum_{k}\left(x_{k}-x_{0}\right)\left(x_{k}-x_{0}\right)^{T} w_{i}} \\ {w_{i}^{T} w_{j}=\delta_{i j} \quad \delta_{i j}=1, \quad i=j} \\ {\delta_{i j}=0, \quad i \neq j}\end{array}
    我们借助拉格朗日乘子法来求解此约束优化问题:
    L=k(xkx0)T(xkx0)i=1dwiTk(xkx0)(xkx0)Twii=1dλi(wiTwi1)L=\sum_{k}\left(x_{k}-x_{0}\right)^{T}\left(x_{k}-x_{0}\right)-\sum_{i=1}^{d} w_{i}^{T} \sum_{k}\left(x_{k}-x_{0}\right)\left(x_{k}-x_{0}\right)^{T} w_{i}-\sum_{i=1}^{d} \lambda_{i}\left(w_{i}^{T} w_{i}-1\right)
    Lx0=2(Ipi=1dwiwiT)k(xkx0)Lwi=2k(xkx0)(xkx0)Twi2λiwi\begin{array}{l}{\frac{\partial L}{\partial x_{0}}=-2\left(I_{p}-\sum_{i=1}^{d} w_{i} w_{i}^{T}\right) \sum_{k}\left(x_{k}-x_{0}\right)} \\ {\frac{\partial L}{\partial w_{i}}=2 \sum_{k}\left(x_{k}-x_{0}\right)\left(x_{k}-x_{0}\right)^{T} w_{i}-2 \lambda_{i} w_{i}}\end{array}
    由两个偏导为0,可以得到:
    x0=kxkNk(xkx0)(xkx0)Twi=λiwi\begin{array}{l}{x_{0}=\sum_{k} \frac{x_{k}}{N}} \\ {\sum_{k}\left(x_{k}-x_{0}\right)\left(x_{k}-x_{0}\right)^{T} w_{i}=\lambda_{i} w_{i}}\end{array}
    因为半正定矩阵的特征值非负,所以,原最小化损失函数可进行转化:
    minXkDsX(xk,X)=k(xkx0)T(xkx0)i=1dwiTk(xkx0)(xkx0)Twi=k(xkx0)T(xkx0)i=1dλiwiTwi=k(xkx0)T(xkx0)i=1dλi\begin{array}{l}{\min _{\underline{X}} \sum_{k} D s_{X}\left(x_{k}, \underline{X}\right)} \\ {=\sum_{k}\left(x_{k}-x_{0}\right)^{T}\left(x_{k}-x_{0}\right)-\sum_{i=1}^{d} w_{i}^{T} \sum_{k}\left(x_{k}-x_{0}\right)\left(x_{k}-x_{0}\right)^{T} w_{i}} \\ {=\sum_{k}\left(x_{k}-x_{0}\right)^{T}\left(x_{k}-x_{0}\right)-\sum_{i=1}^{d} \lambda_{i} w_{i}^{T} w_{i}} \\ {=\sum_{k}\left(x_{k}-x_{0}\right)^{T}\left(x_{k}-x_{0}\right)-\sum_{i=1}^{d} \lambda_{i}}\end{array}
    我们利用矩阵的性质,要想最小化距离平方和,有:
    minXkDsX(xk,X)=k(xkx0)T(xkx0)i=1dλi\min _{\underline{X}} \sum_{k} \operatorname{Ds}_{X}\left(x_{k}, \underline{X}\right)=\sum_{k}\left(x_{k}-x_{0}\right)^{T}\left(x_{k}-x_{0}\right)-\sum_{i=1}^{d} \lambda_{i}
    ΣX=k(xkx0)(xkx0)T\Sigma_{X}=\sum_{k}\left(x_{k}-x_{0}\right)\left(x_{k}-x_{0}\right)^{T}p×pp\times p的矩阵。有性质:
    tr(ΣX)=k(xkx0)T(xkx0)=i=1pλi\operatorname{tr}\left(\Sigma_{X}\right)=\sum_{k}\left(x_{k}-x_{0}\right)^{T}\left(x_{k}-x_{0}\right)=\sum_{i=1}^{p} \lambda_{i}
    则有,
    minXkDsX(xk,X)=i=1pλii=1dλi=i=d+1pλi\min _{\underline{X}} \sum_{k} \operatorname{Ds}_{X}\left(x_{k}, \underline{X}\right)=\sum_{i=1}^{p} \lambda_{i}-\sum_{i=1}^{d} \lambda_{i}=\sum_{i=d+1}^{p} \lambda_{i}
    由此我们可以看到,要得到极小值,我们只要计算XXTXX^T矩阵的前d个最大特征值,是投影后样本具有最小损失的特点。那么此时的WW就是XXTXX^T矩阵前d个最大特征值对应的特征向量。
    不难知道,对于XXTXX^T的特征分解: XXT=UΛUTXX^T = U\Lambda U^T
    这里的U就是前天提到的奇异值分解的U。同理,虽然我们这里没有用到VV,但其实奇异值分解的VV正式XTXX^TX的特征值分解的特征矩阵。
    为了比较XXTXX^T特征分解和XX进行奇异值分解的消耗,写了一段小程序,并使用matlab探查功能进行比较如下:
    在这里插入图片描述

    这个比较事实上没有太大的意义。所用的代码如附录。

    PCA算法简单描述

    假设XX是一个m*n矩阵,表示n个对象的m个特征表示数据,即每一列表示一个对象,每一行表示一个特征。我们希望将特征降为d维,d远小于m。输出结果为YY,一个d*n的矩阵。

    • X=[x1,x2...xn]X=[x_1,x_2...x_n],计算每个对象点的平均值x0=1ni=1nxix_0 = \frac{1}{n}\sum\limits _{i=1}^nx_i

    • Xx0:=[x1x0,x2x0...xnx0]X-x_0 : = [x_1-x_0,x_2-x_0...x_n-x_0]做奇异值分解:Xx0=UΛVTX-x_0 = U\Lambda V^T

    • x0x_0即为新坐标系的原点,UU的前d列即为去中心化后的新的坐标系,不妨记为WW。那么,所有点在新坐标系下的表示为:Y=WT(Xx0)Y=W^T*(X-x_0)。同样地,要将新的投影点yy还原到原坐标系中,可以写为:x0+Wyx_0+W*y

    下面以基于矩阵的视角写出PCA算法的算法流程,输入为矩阵p*N矩阵X,输出为d*N矩阵Y。矩阵的每一列都表示一个对象,每一行都表示对象的一个特征表示。

    在这里插入图片描述

    PCA算例一

    假设小明和小红有身高和体重两个特征(实际操作数据要进行预处理,这里不做),如下表:

    在这里插入图片描述
    那么此时X=[178 165;70 65]X = [178 ~165; 70 ~65],现在试图通过PCA降维,将身高和体重合并为一个特征。走一遍上面的过程,可得:

    Xx0=UΛVTX-x_0 = U\Lambda V^T

    其中,

    那么,有

    那就是说,最后数据可降维为:

    在这里插入图片描述
    这个问题MATLAB计算的小程序在附录。

    PCA算例二

    这是一个对于人脸数据进行降维的例子,人脸数据是我从网上找的。MATLAB源代码见附录。
    选取了2000x1680的数据集进行了测试,选取降维后维数为20,其降维前后的图像(降维后的图像指的是投影点还原到原空间对应的坐标值重构出的图像)如下所示(选取第一个点为代表):
    在这里插入图片描述
    我们使用别人制作的降维工具箱"drtoolbox"重新进行计算并和我的程序结果进行比较。工具箱的使用代码见附录。结果如下:
    在这里插入图片描述
    当然,我们也可以比较我的程序和工具箱程序的误差的大小,比如L2L_2误差。都很简单,暂且不提。

    其他数据降维方法

    其他的数据降维方法还有很多,比如说线性判别分析,拉普拉斯特征映射等等,我这里就简单介绍一下局部线性嵌入。

    当数据具备某些非线性结构,如流形结构时,我们希望降维后的数据仍然保持这些结构。那么就提出了LLE降维算法。LLE(Locally linear embedding):在数据降维后仍然保留原始高维数据的拓扑结构,这种拓扑结构表现为数据点的局部邻接关系。

    此算法我们首先要寻求每个数据点的k个最近邻,然后将当前数据点用k个最近邻线性表出,那么就有相对的权重系数。
    我们希望数据在降维后数据点之间依然能保持这种线性表出的关系,并且在满足另外一些约束条件的前提下,我们很容易求得降维后的数据。
    具体原理和公式网络上有很多人整理得很好,这里不提了。

    下面是LLE算法的算法流程,输入为矩阵p*N矩阵X,输出为d*N矩阵Y。矩阵的每一列都表示一个对象,每一行都表示对象的一个特征表示。
    在这里插入图片描述
    源代码见附录。

    选取了409×698的图像数据集进行了测试,选取降维后维数为2,选取最近邻个数k=12k=12,实验后的部分结果如下:
    在这里插入图片描述
    我们使用别人制作的降维工具箱"drtoolbox"重新进行计算并和我的程序结果进行比较。工具箱的使用代码见附录。

    降维后的部分数据截图如下:
    在这里插入图片描述
    为了比较性能,找个一个别人写的LEE算法,算是网络版本,代码在附录。"网络版"的数据结果和我的版本的结果是一样的。我们开启Matlab的探查功能来比较耗时,结果如图。
    在这里插入图片描述

    函数型数据主成分分析

    Idea的萌生

    前一段时间我在做一个流体力学上的东西(虽然现在已经不做这个方向了),其中比较关键的步骤就是需要用一个带时间变量的多项式公式,来刻画一个物理过程。这个多项式的各个项前面的系数是未知的,由物理规律来决定。我们希望从一些物理实验数据中来通过一些机器学习的手段来学到多项式各个项前面的系数。

    这个问题本质的困难在于,我们不知道那些函数项(基函数)是我们需要的。事实上,只要知道了多项式包含哪些项,是可以通过一些物理原理求得前面的系数的。一个基本的想法就是选足够多的基函数,使得函数空间足够大而包含真值。但是,函数空间太大会带来使用物理原理求系数时的计算困难增大。所以,我们希望能找一个原来大的函数空间的一个子空间,使得用这个子空间,就能够基本刻画原来的物理过程。再用物理原理来求得以子空间基函数为各个项的多项式系数。

    仔细一想,这不正是函数空间的PCA吗?如果把每一个函数看做一个数据点,把各个基函数看做是组成坐标系的坐标轴,那么"函数点"在高维函数空间中的表示,就可以通过类似于主成分分析的技巧,变成在低维函数空间中的表示。只要有了能表示刻画整个物理过程的各个数据点的低维空间,那么刻画物理过程的多项式的项(即低维空间的基函数)也就明确了,剩下的事情也就自然而然了。

    FPCA简介和理论推导

    函数型主成分分析(FPCA,Functional Principal Components Analysis)是传统的PCA的一种推广。考虑我们已经从数据中得到拟合曲线xi(s),sT,i=1,,nx_{i}(s), s \in \mathcal{T}, i=1, \cdots, n,所谓的第一主成分,就是我们希望能找到一个模为1的函数β(s)\beta(s),使得{xi}\{x_i\}β\beta上的投影(L2L_2內积){ξi}\{\xi _i\}的方差达到最大,方差最大其实也就体现{xi}\{x_i\}整体到β\beta的距离达到最小。β\beta一般就叫做权重函数(可以理解为"坐标轴"单位长度量)。

    我们管各个函数到β\beta上的投影叫做观测曲线的主成分得分:
    ξi=Tβ(s)xi(s)ds,i=1,,n\xi_{i}=\int_{\mathcal{T}} \beta(s) x_{i}(s) d s, \quad i=1, \cdots, n故而,求解第一个主成分就变成了求解一个优化问题:
    max1ni=1nξi2=max1ni=1n(Tβ(s)xi(s)ds)2 s.t. β2=Tβ(s)β(s)ds=1\begin{aligned} \max \frac{1}{n} \sum_{i=1}^{n} \xi_{i}^{2} &=\max \frac{1}{n} \sum_{i=1}^{n}\left(\int_{\mathcal{T}} \beta(s) x_{i}(s) d s\right)^{2} \\ \text { s.t. } &\|\beta\|^{2}=\int_{T} \beta(s) \beta(s) d s=1 \end{aligned}求解这个优化问题,我们就得到了第一主成分β1(s)\beta^1(s)
    kk主成分无非就是在满足和前面k1k-1个主成分权重函数垂直的基础上,求解上述优化问题而已,即求解
    max1ni=1nξi2=max1ni=1n(Tβ(s)xi(s)ds)2 s.t. β2=Tβ(s)β(s)ds=1Tβ(s)βl(s)ds=0,l=1,,k1\begin{array}{l}{\max \frac{1}{n} \sum_{i=1}^{n} \xi_{i}^{2}=\max \frac{1}{n} \sum_{i=1}^{n}\left(\int_{\mathcal{T}} \beta(s) x_{i}(s) d s\right)^{2}} \\ {\text { s.t. }\|\beta\|^{2}=\int_{T} \beta(s) \beta(s) d s=1} \\ {\int_{T} \beta(s) \beta^{l}(s) d s=0, l=1, \cdots, k-1}\end{array}
    这个优化问题的解可以表述如下。记协方差函数:
    v(s,t)=1n1i=1n(xi(s)x(s))(xi(t)x(t))v(s, t)=\frac{1}{n-1} \sum_{i=1}^{n}\left(x_{i}(s)-\overline{x}(s)\right)\left(x_{i}(t)-\overline{x}(t)\right)
    那么权重函数满足特征方程:
    Tv(s,t)β(t)dt=λβ(s)\int_{\mathcal{T}} v(s, t) \beta(t) d t=\lambda \beta(s)
    定义积分变换: Vβ(s)=Tv(s,t)β(t)dtV \beta(s)=\int_{\mathcal{T}} v(s, t) \beta(t) d t
    这里的VV称为协方差算子,它将函数β\beta变成一个函数。那么,我们有:
    Vβ(s)=λβ(s)V \beta(s)=\lambda \beta(s)
    我们也类比PCA,使用特征值的累积贡献率来衡量主成分所占比例:
    FVE=i=1Kλi/i=1n1λi\mathrm{FVE}=\sum_{i=1}^{K} \lambda_{i} / \sum_{i=1}^{n-1} \lambda_{i}这里之所以对λ\lambda只累计到nn是因为协方差算子VV的秩为样本数量减一个,则非零特征根的个数最多为n1n-1个。
    由上述已知,我们求解主成分最后归结为求解一个特征值问题。
    求解这个问题,目前比较流行的有三种方法:

    • 对函数进行SVD离散化

    • 对函数进行基函数展开

    • 运用一般性的数值积分方法

    我们最后需要的是特征函数,为了避免插值而带来更大的误差,我选用对基函数进行展开的方法。下面简单介绍一个对函数进行基函数展开的基本思路。
    我们的样本基函数xix_i可以通过基函数展开,如下:
    Xi(s)=k=1KcikΦk(s),i=1,2,,NX_{i}(s)=\sum_{k=1}^{K} c_{i k} \Phi_{k}(s), i=1,2, \ldots, N 我们记
    X=(x1,x2,,xN),Φ=(Φ1,,Φk),C=(cik)N×KX=\left(x_{1}, x_{2}, \ldots, x_{N}\right)^{\prime}, \Phi=\left(\Phi_{1}, \ldots, \Phi_{k}\right)^{\prime}, C=\left(c_{i k}\right)_{N \times K}
    那么样本函数就可以写为等价的矩阵形式X=CΦX=C \Phi。那么协方差函数就可以写为(假设已经标准化):
    v(s,t)=1n1Φ(s)CCΦ(t)v(s, t)=\frac{1}{n-1} \Phi^{\prime}(s) C^{\prime} C \Phi(t)
    定义K阶对称矩阵 W=ΦΦW=\int \Phi \Phi^{\prime}
    当选择正交基的时候,比如说正交傅里叶基,这就是一个单位矩阵。关于这个基如何选取,我们后面还会详谈。
    同样地,将特征函数进行展开:
    β(s)=k=1KbkΦk(s)=Φ(s)b\beta(s)=\sum_{k=1}^{K} b_{k} \Phi_{k}(s)=\Phi^{\prime}(s) b
    将其代入 Tv(s,t)β(t)dt=λβ(s)\int_{\mathcal{T}} v(s, t) \beta(t) d t =\lambda \beta(s)
    就可以得到(N=n1N=n-1):
    1NΦ(s)CCWb=λΦ(s)b\frac{1}{N} \Phi^{\prime}(s) C^{\prime} C W b=\lambda \Phi^{\prime}(s) b
    进一步能得到1NCCWb=λb\frac{1}{N} C^{\prime} C W b=\lambda b,由特征向量正交和单位长度的约束要求,有bkWbk=1,bkWbm=0,kmb_{k}^{\prime} W b_{k}=1, b_{k}^{\prime} W b_{m}=0,k \neq m
    WW做cholesky分解,可得W=LLW=LL'
    定义u=Lbu=L'b,那么上述问题就变成了对称矩阵的代数特征值问题:
    1NLCCLu=λu\frac{1}{N} L' C^{\prime} C L u=\lambda u
    据此可以求得uu,进而求得bb,最后求得特征函数β\beta

    常用的基函数有傅里叶基函数和B样条基函数,傅里叶基函数适用于周期性函数数据,B样条基函数适用于非周期函数数据,当然,也可以用多项式基函数。
    B样条基函数的递归定义为:
    Bj,0(x)={1,tjx<tj+10,elseBi,k(x)=xtiti+ktiBi,k1(x)+ti+k+1xti+k+1ti+1Bi+1,k1(x),k>0\begin{array}{c}{B_{j, 0}(x)=\left\{\begin{array}{l}{1, t_{j} \leq x<t_{j+1}} \\ {0, \text {else}}\end{array}\right.} \\ {B_{i, k}(x)=\frac{x-t_{i}}{t_{i+k}-t_{i}} B_{i, k-1}(x)+\frac{t_{i+k+1}-x}{t_{i+k+1}-t_{i+1}} B_{i+1, k-1}(x), k>0}\end{array}
    附录中有一段简单的以多项式为基的MATLAB代码。

    FPCA和PCA的区别和联系

    如上所述,可以看出,如果所选的基函数是正交的,本质上和PCA的以拟合系数为坐标点的函数空间PCA推广是实际上是一样的。若基函数不是正交的,无非就是在此基础上对要求特征值的矩阵得多乘一个W=ΦΦW=\int \Phi \Phi^{\prime},再求特征向量,以及进行WW意义下对特征向量进行单位化而已(不单位化也没事,只不过权重函数β\beta不再是模长为1的而已,WW意义下的单位话也就意味着让新的基函数模长为1)。这个也非常容易理解,因为在从函数的元(primal)表示左乘一个质量矩阵就到了到它的对偶(dual)表示,而在基函数不正交的情况下,我们应该在对偶空间中再进行它的主成分分析降维,即各个函数的向量表示应该为这个函数和各个基函数的內积。同理,在对偶框架下得到的新的基函数的向量表示也是在对偶空间下的,应该左乘一个质量矩阵才能回到元空间中去。

    基于FPCA的模型约化

    Onsager原理简介

    Onsager基本原理是基于物理规律的一个原理,利用它不难得到,如果刻画物理过程的模型方程有哪些项知道了,也就是基函数知道了,那么我们可以通过这个原理求得各个项前面的系数。
    定义势能函数(自由能): A(a)A(a) 定义能量耗散函数:
    Φ(a˙,a)=12i,jζij(a)a˙ia˙j\Phi(\dot{a}, a)=\frac{1}{2} \sum_{i, j} \zeta_{i j}(a) \dot{a}_{i} \dot{a}_{j}
    那么系统随时间演化由最小化以下函数得到:
    R(a˙,a)=Φ(a˙,a)+iAaia˙iR(\dot{a}, a)=\Phi(\dot{a}, a)+\sum_{i} \frac{\partial A}{\partial a_{i}} \dot{a}_{i}
    最小化RR,可以得到:
    Φa˙i+Aai=0 or jζij(a)a˙j=Aai\frac{\partial \Phi}{\partial \dot{a}_{i}}+\frac{\partial A}{\partial a_{i}}=0 \quad \text { or } \quad \sum_{j} \zeta_{i j}(a) \dot{a}_{j}=-\frac{\partial A}{\partial a_{i}}
    这就是我们要求解的ODE系统。

    简单例子:斜板液滴滑动

    问题描述

    考虑一个液滴在斜板上从静止开始下滑,如图。
    在这里插入图片描述
    从正面或者侧面拍摄到的图案大概如图。
    在这里插入图片描述
    我们现在想要刻画这个液滴的状态,即在每一时刻液滴的俯视形状以及侧视高度。
    我们可以用一个方程来描述这个过程:
    h(x,y,t)=H(x,t)[1(yY(x,t))2]h(x, y, t)=H(x, t)\left[1-\left(\frac{y}{Y(x, t)}\right)^{2}\right]
    其中xx为平行平板沿着液滴运动的方向,yy为平行平板垂直于液滴运动的方向,tt为时刻,hh为垂直于平板距离平板的一个高度。这里面的HYH、Y是两个函数,分别刻画了俯视的形状和侧视的形状。事实上,取h=0h=0,可以得到y=Y(x,t)y=Y(x,t)描述了俯视图(垂直于板)的形状(一半),取y=0y=0,得到h=H(x,t)h=H(x,t),体现的是侧视图。再者,若给定了xx值,高度随着yy是呈现出抛物的变化。因此,这个公式看起来不无道理。
    接下来,我们对H,YH,Y做一个简单的假定:
    H(x,t)=(xa1(t))(a2(t)x)(a3(t)+a4(t)x)Y(x,t)=(xa1(t))12(a2(t)x)12(a5(t)+a6(t)x)\begin{array}{c}{H(x, t)=\left(x-a_{1}(t)\right)\left(a_{2}(t)-x\right)\left(a_{3}(t)+a_{4}(t) x\right)} \\ {Y(x, t)=\left(x-a_{1}(t)\right)^{\frac{1}{2}}\left(a_{2}(t)-x\right)^{\frac{1}{2}}\left(a_{5}(t)+a_{6}(t) x\right)}\end{array}
    容易想到,这里的a1(t),a2(t)a_1(t),a_2(t)表示的是液滴的前后端点(采用欧拉坐标系),因为H,YH,Y在两端点处的值为零。

    原理的应用

    我们希望能通过上面提到的Onsager原理来确定这里的系数aia_i
    固定时刻的液滴体积:
    Ω=a1a2dxYYdyh(x,y,t)\Omega=\int_{a_{1}}^{a_{2}} d x \int_{-Y}^{Y} d y h(x, y, t)
    因为体积是守恒量,所以问题的自由度个数就变成了5。 势能函数定义为:
    A(a)=a1a2dxYYdy[12γθe2+12γ[(xh)2+(yh)2]+12ρgh2sinαρgxhcosα]\begin{aligned} A(a)=& \int_{a_{1}}^{a_{2}} d x \int_{-Y}^{Y} d y\left[\frac{1}{2} \gamma \theta_{e}^{2}+\frac{1}{2} \gamma\left[\left(\partial_{x} h\right)^{2}+\left(\partial_{y} h\right)^{2}\right]\right.\\ &+\frac{1}{2} \rho g h^{2} \sin \alpha-\rho g x h \cos \alpha ] \end{aligned}
    这里的γ\gamma表示液滴的表面张力,ρ\rho表示密度,θe\theta_e是平衡态下的接触角大小,gg是重力加速度,α\alpha是前面提到的斜面角。我也不知道势能函数为什么能写成这样,需要一些物理的分析。
    可以把hh的表达式代入到这个势能函数的表达式中。我们还需要知道能量耗散函数Φ\Phi。由滑润近似,能量耗散函数可以写成关于速度的变量:
    Φ[vx,vy]=12a1a2dxYYdy3ηh(vx2+vy2)\Phi\left[v_{x}, v_{y}\right]=\frac{1}{2} \int_{a_{1}}^{a_{2}} d x \int_{-Y}^{Y} d y \frac{3 \eta}{h}\left(v_{x}^{2}+v_{y}^{2}\right)
    这里的vx,vyv_x,v_y表示两个方向上的速度,η\eta表示流体的粘性。但是我们想要的耗散函数是关于a˙\dot a
    的,所以要想办法替换掉速度。 由体积守恒,我们有:
    h˙=x(vxh)y(vyh)\dot{h}=-\partial_{x}\left(v_{x} h\right)-\partial_{y}\left(v_{y} h\right)
    hh的表达式代入上式,可得:
    (1y2Y2)(H˙+x(vxH)+Hyvy)+2HyY3(yY˙+yvxxYYvy)=0\begin{array}{c}{\left(1-\frac{y^{2}}{Y^{2}}\right)\left(\dot{H}+\partial_{x}\left(v_{x} H\right)+H \partial_{y} v_{y}\right)} \\ {+\frac{2 H y}{Y^{3}}\left(y \dot{Y}+y v_{x} \partial_{x} Y-Y v_{y}\right)=0}\end{array}
    这个约束满足的一个充分条件是:
    H˙+x(vxH)+Hyvy=0yY˙+yvxxYYvy=0\begin{array}{l}{\dot{H}+\partial_{x}\left(v_{x} H\right)+H \partial_{y} v_{y}=0} \\ {y \dot{Y}+y v_{x} \partial_{x} Y-Y v_{y}=0}\end{array}
    一个如下所示的速度场能够满足这样的条件:
    vx(x,y,t)=V(x,t),vy(x,y,t)=W(x,t)yv_{x}(x, y, t)=V(x, t), \quad v_{y}(x, y, t)=W(x, t) y
    其中,V,WV,W的表达为:
    V(x,t)=1HYa1x(H˙Y+HY˙)dxW=1Y(Y˙+VxY)\begin{aligned} V(x, t) &=-\frac{1}{H Y} \int_{a_{1}}^{x}(\dot{H} Y+H \dot{Y}) d x \\ W &=\frac{1}{Y}\left(\dot{Y}+V \partial_{x} Y\right) \end{aligned}
    那么,我们得到的能量耗散函数其实是:
    Φ[a˙,a]=12a1a2dxYYdy3ηh(V2+y2W2)\Phi\left[\dot a, a\right]=\frac{1}{2} \int_{a_{1}}^{a_{2}} d x \int_{-Y}^{Y} d y \frac{3 \eta}{h}\left(V^{2}+{y}^{2}W^2\right)
    我们把aa看成常量,由于H˙,Y˙\dot H,\dot Ya˙\dot a的线性组合,意味着V,WV,W也是,那么Φ\Phi就是a˙\dot a的二次函数,不妨记为:
    Φ(a˙)=12i,jζija˙ia˙j\Phi(\dot{a})=\frac{1}{2} \sum_{i, j} \zeta_{i j} \dot{a}_{i} \dot{a}_{j}
    这里的ξij\xi_{ij}aa的函数。
    这下有了势能函数和能量耗散函数,我们可以得到关于aia_i的发展方程为:
    j=16ζija˙j+Aai=0\sum_{j=1}^{6} \zeta_{i j} \dot{a}_{j}+\frac{\partial A}{\partial a_{i}}=0
    求解之,可得aa

    算法步骤

    总结一下上述的计算过程,就是:

    • 能量耗散函数:
      Φ[a˙,a]=12a1a2dxYYdy3ηh(V2+y2W2)\Phi\left[\dot a, a\right]=\frac{1}{2} \int_{a_{1}}^{a_{2}} d x \int_{-Y}^{Y} d y \frac{3 \eta}{h}\left(V^{2}+{y}^{2}W^2\right)
      其中,
      V(x,t)=1HYa1x(H˙Y+HY˙)dxW=1Y(Y˙+VxY)\begin{aligned} V(x, t) &=-\frac{1}{H Y} \int_{a_{1}}^{x}(\dot{H} Y+H \dot{Y}) d x \\ W &=\frac{1}{Y}\left(\dot{Y}+V \partial_{x} Y\right) \end{aligned}
      h(x,y,t)=H(x,t)[1(yY(x,t))2]h(x, y, t)=H(x, t)\left[1-\left(\frac{y}{Y(x, t)}\right)^{2}\right]
      H(x,t)=(xa1(t))(a2(t)x)(a3(t)+a4(t)x)Y(x,t)=(xa1(t))12(a2(t)x)12(a5(t)+a6(t)x)\begin{array}{c}{H(x, t)=\left(x-a_{1}(t)\right)\left(a_{2}(t)-x\right)\left(a_{3}(t)+a_{4}(t) x\right)} \\ {Y(x, t)=\left(x-a_{1}(t)\right)^{\frac{1}{2}}\left(a_{2}(t)-x\right)^{\frac{1}{2}}\left(a_{5}(t)+a_{6}(t) x\right)}\end{array}

    • 由此,我们计算出Φ(a˙)\Phi(\dot a)表达式,并提取前面的线性组合的系数:
      Φ(a˙)=12i,jζija˙ia˙j\Phi(\dot{a})=\frac{1}{2} \sum_{i, j} \zeta_{i j} \dot{a}_{i} \dot{a}_{j}

    • 势能函数:
      A(a)=a1a2dxYYdy[12γθe2+12γ[(xh)2+(yh)2]+12ρgh2sinαρgxhcosα]\begin{aligned} A(a)=& \int_{a_{1}}^{a_{2}} d x \int_{-Y}^{Y} d y\left[\frac{1}{2} \gamma \theta_{e}^{2}+\frac{1}{2} \gamma\left[\left(\partial_{x} h\right)^{2}+\left(\partial_{y} h\right)^{2}\right]\right.\\ &+\frac{1}{2} \rho g h^{2} \sin \alpha-\rho g x h \cos \alpha ] \end{aligned}

    • 求解ODE方程组(数值解),得出aa
      j=16ζija˙j+Aai=0\sum_{j=1}^{6} \zeta_{i j} \dot{a}_{j}+\frac{\partial A}{\partial a_{i}}=0

    数值实验

    所用的参数如下:
    η=104cP,ρ=964kgm3\eta=104 \mathrm{cP}, \rho=964 \mathrm{kg} \mathrm{m}^{-3}
    γ=20.9mNm1,\gamma=20.9 \mathrm{mNm}^{-1}, θe=53\theta_{e}=53^{\circ}
    Ω=6.3mm3\Omega = 6.3 \mathrm{mm}^{3} α=15,25,45\alpha=15^{\circ},25^\circ,45^\circ

    下面有一些数值结果如图。
    在这里插入图片描述
    所用的程序比较冗长,就不往本文后面贴了。

    FPCA在液滴下滑问题的应用

    这只是我的一个想法,目前有很多问题都没有明确。由于时间关系,我这里也不会展开细述这一部分内容。基本的做法可以分成以下几个步骤:

    • 收集数据:除了网络上搜到的三个物理实验视频和论文中的一些截图之外,我没有找到更多的数据,数据严重不足。和文章作者联系,也未要到数据。

    • 图像处理:对收集到的视频,按帧提取图像,对每个图像进行去噪,二值化,归一化,提取边缘的坐标位置。

    • FPCA降维:对于提取到的数据,选用适当的基函数,做小二乘意义下的拟合,得到拟合系数。这一组组拟合系数,就是我们做FPCA降维的数据。做FPCA,得到子函数空间。

    • Onsager原理确定系数:在子函数空间中,使用Onsager基本原理,得到液滴下滑物理过程的表达式系数。

    数据不够怎么办?有两个基本的想法。一个是利用同一组参数(如斜板角度)下不同时刻的数据(一个视频),来降维生成这组参数下的随时间变化的物理过程表达过程。另一个是查找更多的数据,哪怕利用上别人文章中的图片,堆砌所有的数据,寻求刻画这个物理过程的"真"表达,找到物理上的"真"规律。

    收集到的原始数据如图所示。
    在这里插入图片描述
    处理后的数据如图所示。
    在这里插入图片描述在这里插入图片描述
    其中用到的一些代码见附录。根据这个问题的特殊性,有一个新的想法就是Robust
    PCA和流行学习能不能推广到FPCA上?这也是一个有趣的问题。其实我们还是不太清楚这个问题中数据的分布。

    参考文献

    [1] Rio E , Daerr A , Andreotti B , et al. Boundary Conditions in the
    Vicinity of a Dynamic Contact Line: Experimental Investigation of
    Viscous Drops Sliding Down an Inclined Plane[J]. Physical Review
    Letters, 2005, 94(2):024503.

    [2] Rudy S H , Brunton S L , Proctor J L , et al. Data-driven
    discovery of partial differential equations[J]. Science Advances,
    2017, 3(4):e1602614.

    [3] Brunton S L , Proctor J L , Kutz J N . Discovering governing
    equations from data by sparse identification of nonlinear dynamical
    systems[J]. Proceedings of the National Academy of Sciences,
    2016:201517384.

    [4] Xu X , Di Y , Doi M . Variational method for liquids moving on a
    substrate[J]. Physics of Fluids, 2016, 28(8):087101.

    [5] 胡宇. 函数型数据分析方法研究及其应用[D]. 东北师范大学, 2011.

    [6] 陈宜治. 函数型数据分析若干方法及应用[D]. 浙江工商大学, 2011.

    [7] 沈关友. 基于函数型数据主成分分析的银行股票数据预测[D].

    [8] 吴刚, 胡新荣. 基于函数型主成分分析的织物状态研究[J].
    科技创业月刊, 2017(12).

    [9] 李敏. 基于函数型主成分分析方法的用水量数据分析[J].
    合肥学院学报(综合版), 2014(4):21-25.
    ear dynamical
    systems[J]. Proceedings of the National Academy of Sciences,
    2016:201517384.

    [4] Xu X , Di Y , Doi M . Variational method for liquids moving on a
    substrate[J]. Physics of Fluids, 2016, 28(8):087101.

    [5] 胡宇. 函数型数据分析方法研究及其应用[D]. 东北师范大学, 2011.

    [6] 陈宜治. 函数型数据分析若干方法及应用[D]. 浙江工商大学, 2011.

    [7] 沈关友. 基于函数型数据主成分分析的银行股票数据预测[D].

    [8] 吴刚, 胡新荣. 基于函数型主成分分析的织物状态研究[J].
    科技创业月刊, 2017(12).

    [9] 李敏. 基于函数型主成分分析方法的用水量数据分析[J].
    合肥学院学报(综合版), 2014(4):21-25.

    附录

    {a}
    
    clc
    clear
    p = 10000;
    N = 10;
    d = 10;
    AR = randn(p,N);
    x0 = mean(AR,2); % 计算样本均值
    AR_shift = AR - repmat(x0,1,N); % 中心平移每个训练样本
    Sigma = AR_shift*AR_shift';
    [W,D] = eigs(Sigma,d); % 前d个特征向量作为wi
    Lambda = diag(D); % 提取特征值
    %%
    [U,S,V] = svd(AR_shift);
    
    
    {b}
    \begin{lstlisting}%[language=MATLAB,frame=shadowbox]
    clc
    clear
    X = [178 165; 70 65];
    x0 = sum(X,2)/2;
    X_x0 = X - [x0 x0];
    [U,S,V] = svd(X_x0);
    W = U(:,1);
    Y = W'*(X_x0);
    rebuild_X = [x0 x0]+W*Y;
    
    {c}
    \begin{lstlisting}%[language=MATLAB,frame=shadowbox]
    %%
    %编写程序,实现PCA算法
    %对图像进行降维实验,并显示降维重建后的图像
    %运行已有程序,和自己的对比
    %实验报告(伪代码(或流程图)、源代码、实验结果及分析)
    %%
    %编写程序,实现PCA和LEE算法
    %对图像进行降维实验,并显示降维重建后的图像
    %运行已有程序,和自己的对比
    %实验报告(伪代码(或流程图)、源代码、实验结果及分析)
    clc
    clear
    addpath(genpath(pwd));%将子孙文件添加到工作目录
    load AR %导入数据
    d = 20;
    %%
    AR = double(AR);%双进度化
    [p ,N] = size(AR);%特征维度和对象数目
    % for i = 1:N
    % image = AR(:,i);
    % image = reshape(image,50,40);
    % imshow(mat2gray(image));%对原矩阵归一化
    % end
    x0 = mean(AR,2);%计算样本均值
    AR_shift = AR - repmat(x0,1,N);%中心平移每个训练样本
    %计算协方差矩阵的(n-1)倍,不用内置函数cov,提高代码的重用率和运行速度
    Sigma = AR_shift*AR_shift';
    %Sigma = cov(AR')*(N-1);
    [W,D] = eigs(Sigma,d);%前d个特征向量作为wi
    Lambda = diag(D);%提取特征值
    %%
    close all;
    k = 1;
    Y = W'*AR_shift(:,k);%第k个图像的输出表示
    X_rebuid = W*Y + x0;%第k个图像的重建还原
    
    image = AR(:,k);
    image = reshape(image,50,40);
    imshow(mat2gray(image));%对原矩阵归一化
    
    figure;
    image_re = X_rebuid;
    image_re = reshape(image_re,50,40);
    imshow(mat2gray(image_re));%对原矩阵归一化
    
    
    
    {d}
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    %%
    %编写程序,实现PCA和LLE算法
    %对图像进行降维实验,并显示降维重建后的图像
    %运行已有程序,和自己的对比
    %实验报告(伪代码(或流程图)、源代码、实验结果及分析)
    %%
    %编写程序,实现PCA和LEE算法
    %对图像进行降维实验,并显示降维重建后的图像
    %运行已有程序,和自己的对比
    %实验报告(伪代码(或流程图)、源代码、实验结果及分析)
    clc
    clear
    addpath(genpath(pwd));%将子孙文件添加到工作目录
    load AR %导入数据
    d = 20;
    
    %%
    AR = double(AR);%双进度化
    [p ,N] = size(AR);%特征维度和对象数目
    % for i = 1:N
    % image = AR(:,i);
    % image = reshape(image,50,40);
    % imshow(mat2gray(image));%对原矩阵归一化
    % end
    x0 = mean(AR,2);%计算样本均值
    AR_shift = AR - repmat(x0,1,N);%中心平移每个训练样本
    %计算协方差矩阵的(n-1)倍,不用内置函数cov,提高代码的重用率和运行速度
    Sigma = AR_shift*AR_shift';
    %Sigma = cov(AR')*(N-1);
    [W,D] = eigs(Sigma,d);%前d个特征向量作为wi
    Lambda = diag(D);%提取特征值
    %%
    close all;
    k = 1;
    Y = W'*AR_shift(:,k);%第k个图像的输出表示
    X_rebuid = W*Y + x0;%第k个图像的重建还原
    
    image = AR(:,k);
    image = reshape(image,50,40);
    imshow(mat2gray(image));%对原矩阵归一化
    
    figure;
    image_re = X_rebuid;
    image_re = reshape(image_re,50,40);
    imshow(mat2gray(image_re));%对原矩阵归一化
    
    
    
    {e}
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    %%
    %编写程序,实现PCA和LEE算法
    %对图像进行降维实验,并显示降维重建后的图像
    %运行已有程序,和自己的对比
    %实验报告(伪代码(或流程图)、源代码、实验结果及分析)
    
    %% 预处理和数据输入
    clc
    clear
    addpath(genpath(pwd));%将子孙文件添加到工作目录
    load face_images; %导入数据
    data = images;
    %data = data(:,1:50);
    %% 初始化参数
    d = 2;
    len = 64;
    wid = 64;
    k = 12;
    %%
    [p ,N] = size(data);%特征维度和对象数目
    [IDX,~] = knnsearch(data',data','K',k+1);
    IDX = IDX(:,2:end);
    W = zeros(N);
    for i = 1:N
        xk = data(:,i);
        index = IDX(i,:);
        Qk_temp = repmat(xk,1,k) - data(:,index);
        Qk = Qk_temp'*Qk_temp;
        wk_temp = Qk\ones(k,1);
        wk = wk_temp/sum(wk_temp);
        W(i,index) = wk;
    end
    W = W';
    I = eye(N);
    M = (I-W)*(I-W)';
    [P,L] = eigs(M,d+1,0);
    P = P(:,2:end);
    Y = (P*sqrt(N))';
    
    
    
    
    
    
    
    {f}
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    %% 使用工具箱进行进行降维来和我的实验结果进行比较
    clc
    clear
    close all
    
    method = 'LLE';%可选LLE或者PCA
    addpath(genpath(pwd));
    
    % 产生测试数据
    %[X, labels] = generate_data('helix', 2000);
    if strcmp(method,'PCA')
        load AR %导入数据
        [p,N] = size(AR);
        X = double(AR);%导入数据
    else
        load face_images %导入数据
        [p,N] = size(images);
        X = double(images);%导入数据
    end
    
    
    % 估计本质维数
    %no_dims = round(intrinsic_dim(X, 'MLE'));
    %disp(['MLE estimate of intrinsic dimensionality: ' num2str(no_dims)]);
    d = 2;
    k = 12;
    % PCA降维或LLE降维
    [mappedX, mapping] = compute_mapping(X', method,d);
    Y = mappedX';
    
    if strcmp(method,'PCA')
        x0 = (mapping.mean)';
        W = (mapping.M);
        AR_shift = X - repmat(x0,1,N);
    
        %%
        close all;
        k = 1;
        y = Y(:,k);
        X_rebuid = W*y + x0;%第k个图像的重建还原
    
        image = AR(:,k);
        image = reshape(image,50,40);
        imshow(mat2gray(image));%对原矩阵归一化
    
        figure;
        image_re = X_rebuid;
        image_re = reshape(image_re,50,40);
        imshow(mat2gray(image_re));%对原矩阵归一化
    
    end
    
    
    
    {g}
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    % LLE ALGORITHM (using K nearest neighbors)
    % [Y] = lle(X,K,dmax)
    % X    :data as D x N matrix (D = dimensionality, N = #points)
    % K    :number of neighbors
    % dmax :max embedding dimensionality
    % Y    :embedding as dmax x N matrix
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %function [Y] = lle(X,K,d)
    addpath(genpath(pwd));%将子孙文件添加到工作目录
    load face_images; %导入数据
    data = images;
    X = data;
    K = 12;
    d = 2;
    %%
    [D,N] = size(X);
    fprintf(1,'LLE running on %d points in %d dimensions\n',N,D);
    
    %% Step1: compute pairwise distances & find neighbour
    fprintf(1,'-->Finding %d nearest neighbours.\n',K);
    
    X2 = sum(X.^2,1);
    distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;
    [sorted,index] = sort(distance);
    neighborhood = index(2:(1+K),:);
    
    % Step2: solve for recinstruction weights
    fprintf(1,'-->Solving for reconstruction weights.\n');
    if(K>D)
      fprintf(1,'   [note: K>D; regularization will be used]\n');
      tol=1e-3; % regularlizer in case constrained fits are ill conditioned
    else
      tol=0;
    end
    
    W = zeros(K,N);
    for ii=1:N
       z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); % shift ith pt to origin
       C = z'*z;                                        % local covariance
       C = C + eye(K,K)*tol*trace(C);                   % regularlization (K>D)
       W(:,ii) = C\ones(K,1);                           % solve Cw=1
       W(:,ii) = W(:,ii)/sum(W(:,ii));                  % enforce sum(w)=1
    end;
    
    % Step 3: compute embedding from eigenvects of cost matrix M=(I-W)'(I-W)
    fprintf(1,'-->Computing embedding.\n');
    % M=eye(N,N); % use a sparse matrix with storage for 4KN nonzero elements
    M = sparse(1:N,1:N,ones(1,N),N,N,4*K*N);
    for ii=1:N
       w = W(:,ii);
       jj = neighborhood(:,ii);
       M(ii,jj) = M(ii,jj) - w';
       M(jj,ii) = M(jj,ii) - w;
       M(jj,jj) = M(jj,jj) + w*w';
    end;
    
    % calculation of embedding
    options.disp = 0;
    options.isreal = 1;
    options.issym = 1;
    [Y,eigenvals] = eigs(M,d+1,0,options);
    Y = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0
    fprintf(1,'Done.\n');
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % other possible regularizers for K>D
    %   C = C + tol*diag(diag(C));                       % regularlization
    %   C = C + eye(K,K)*tol*trace(C)*K;                 % regularlization
    
    
    
    
    
    {h}
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    clc
    clear
    close all
    format short
    digits(5)
    %% 设定维数
    p = 8;
    d = 2;
    %% 读入数据,保存
    D = dir('../data/*.csv');
    D = struct2cell(D);%结构体不好调用,转成元胞数组来使用
    csvs_name = D(1,:);
    N = 0;%这个表示数据的个数
    %csvs_name = csvs_name(5:6);
    for csv_name = csvs_name %cpicName = 'a001.bmp'
        N = N+1;
        csvName = csv_name{1};
        Data{N} = load(['../data/' csvName]);
    end
    %clear csv_name csvN csvName csvs_name D
    %% 定义原空间基,生成模型函数,由模型函数通过拟合得到点的坐标表示
    %p = 15;%原空间维数
    %%根据p自动生成模型函数
    eval_poly = 'a(1)';
    for i = 2:p
        eval_poly = [eval_poly '+a(' num2str(i) ').*x.^' num2str(i-1)];
    end
    eval_poly = ['modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(' eval_poly ');'];
    eval(eval_poly);
    
    %modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*sin(x)+a(3).*cos(x)+a(4).*1./x+a(5).*x+a(6).*x.^2);
    %   modelfun = @(a,x)(x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*x+a(3).*x.^2+...
    %       a(4).*x.^3))%+a(5).*x.^4+a(6)*x.^5))+...
    %       a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
    %         modelfun = @(a,x)((a(1)+a(2).*x.^2+a(3).*x.^2+...
    %         a(4).*x.^3+a(5).*x.^4+a(6)*x.^5+...
    %         a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
    
    %%做一波拟合,求出系数,进行拟合,拟合结果为A
    %%A中存的是数据点的坐标表示
    for k = 1:N
        data = Data{k};
        data(:,2) = abs(data(:,2));%取绝对值
        x = data(:,1);
        y = data(:,2);
        %y = y./(x.*(1-x));%%%%%%%%如果将x(1-x)项挪到左边,拟合会出问题
        a0 = ones(1,p);
        a = nlinfit(x,y,modelfun,a0);
        A(:,k) = a';
    end
    %clearvars -except A modelfun x y
    
    %% 开始奇异值分解降维,得到新空间的原点和坐标轴,坐标表示
    %%由奇异值分解做主成分分析,A的每一列为原空间下的系数,U的每一列为主成分
    [p,N] = size(A);
    %d = 5;%新空间维数
    x0 = mean(A,2);%计算样本均值
    A_shift = A - repmat(x0,1,N);%中心平移每个训练样本,A_shift为原数据点的平移表示
    
    %%%%对于A_shift的补充处理
    syms x;
    base_vec_T = x.^(0.5).*(1-x).^(0.5).*x.^[0:p-1];
    base_vec = conj((base_vec_T))';
    W = int(base_vec*base_vec_T,[0,1]);
    L_prime = chol(W);
    L_prime_A_shift = L_prime*A_shift;
    
    %%%%
    [U,S,V] = svd(L_prime_A_shift);
    B = L_prime\U;
    if size(S,2) == 1
        S_diag = S(1);
    else
        S_diag = diag(S);
    end
    R = cumsum(S_diag)./sum(S_diag);
    Coord_new = B(:,1:d);
    
    
    %% 新空间的坐标轴和原点的函数表示
    %%生成新的空间中的基函数,找到d个新空间基函数及原点坐标
    %%Coord_new x0中存的是新坐标系,fi存的是新的函数空间的原点和坐标轴
    for k = 1:d;%取出第k个基函数
        syms x;
        c = B(:,k);
        exp = modelfun(c,x);
        exp = simplify(exp);%第k个基函数的表达式
        f{k} = exp;
        digits(2);
        latex(vpa(exp,2));%写成latex表达式黏贴到解释器中
    end
    
    f0 = modelfun(x0,x);
    f0 = simplify(f0);
    latex(vpa(f0,2));
    % for i = 1:d
    %     f_fun{i} = matlabFunction(f{i});
    % end
    % f0_fun = matlabFunction(f0);
    
    %% 求原坐标系中的投影后在原空间中的位置
    
    
    
    %A_new_coord = Coord_new'*A_shift;%A_new_coord表示原数据点在新空间中的表示
    
    syms x;
    for k = 1:N
        a = A_shift(:,k);
        pk_f = modelfun(a,x);
        for j = 1:d
            A_new_coord(j,k) = int(pk_f.*f{j},0,1);
        end
    end
    
    %%函数运算
    for k = 1:N
        A_shadow_f = f0;
        for i = 1:d
            A_shadow_f = A_shadow_f+f{i}*A_new_coord(i,k);
        end
        A_shadow_f = simplify(A_shadow_f);
        A_shadow_fs{k} = A_shadow_f;
        A_shadow_fs_fun{k} = matlabFunction(A_shadow_f);
    end
    
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % %%坐标运算,可以注释掉
    % A_shadow_coord = Coord_new*A_new_coord+repmat(x0,1,N);%A_shadow_coord表示投影坐标点在原坐标系下的表示
    % syms x;
    % for k = 1:N
    %     a = A_shadow_coord(:,k);
    %     A_shadow_coord_f = modelfun(a,x);
    %     A_shadow_coord_fs{k} = simplify(A_shadow_coord_f);
    % end
    %
    % %%不考虑舍入误差的情况下,可以证明坐标运算和函数运算的得到的函数表达式是一样的,验证代码如下
    % k=5;
    % diff = A_shadow_coord_fs{k}-A_shadow_fs{k}
    % diff = simplify(diff)
    % digits(2);
    % latex(vpa(diff,2))
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %% 画图和求误差
    Errors = [];
    NN = 1000;
    h = 1/NN;
    xx = linspace(0,1,NN+1);
    close all;
    for k = 1:N
        %%原数据点
        data = Data{k};
        data(:,2) = abs(data(:,2));%取绝对值
        x = data(:,1);
        y = data(:,2);
        %%拟合曲线
        a = A(:,k);
        data_fit_f = @(x) modelfun(a,x);
        %%投影函数
        A_shadow_f_fun = A_shadow_fs_fun{k};
        %%原点函数
        f0_fun = matlabFunction(f0);
    
        %%开始画图
        figure(k);
        scatter(x,y,'.','MarkerEdgeColor','b',...
                  'MarkerFaceColor','c',...
                  'LineWidth',0.5);
        hold on;
    
        y_fit = data_fit_f(xx);
        y_shadow = A_shadow_f_fun(xx);
        y_f0 = f0_fun(xx);
        plot(xx,y_fit,'g');
        plot(xx,y_shadow,'black');
        plot(xx,y_f0,'r');
        legend('数据','拟合','投影','原点');
    
        %%求拟合曲线和投影函数之间的L2误差
        %Errors(end+1) = norm((y_fit-y_shadow),2);
        Errors(end+1) = sum(abs(y_fit-y_shadow)*h);
    end
    
    Errors_mean = mean(Errors)
    %% test
    % close all;
    % plot(xx,xx);
    % hold on;
    % plot(xx,xx.^9)
    
    
    {i}
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    %图片提取数据与清洗
    function I_fresh = denoise( I )
    % 去除噪声,填充空洞
    [L,num] = bwlabel(~I,8);
    stats = regionprops(L);
    stats_cell = struct2cell(stats);
    [~,maxind] = max([stats_cell{1,:}]);
    Imp = ismember(L, maxind);
    I_fresh = imfill(Imp,'holes');
    I_fresh = ~I_fresh;
    end
    
    function [ coord,L ] = edge_extract( I )
    % 提取边缘
    [B,L] = bwboundaries(I,'noholes');
    %imshow(I)
    % figure;
    % imshow(L)
    coord = B{1};
    L = bwperim(I);
    end
    
    function [Ia,Ib] = image_process( pic_name,plotflag)
    % 图像分割,二值化处理
    %pic_name = '1\65.png';
    A = imread(pic_name);%读取到一张图片
    [low_num,col_num,~] = size(A);
    Aa = imcrop(A,[1,1,270-1,low_num-1]);%分割点定在了270
    Ab = imcrop(A,[270,1,col_num-270,low_num-1]);
    %thresh = graythresh(A);%自动确定二值化阈值
    %L = bwperim(A);imshow(L)
    Ia = im2bw(Aa,58/255); %对图像二值化处理
    Ib = im2bw(Ab,58/255); %对图像二值化处理
    
    
    if plotflag == 1
    %if 0
        figure(100);
        subplot(2,2,1);imshow(Aa);
        subplot(2,2,2);imshow(Ab);
        subplot(2,2,3);imshow(Ia);
        subplot(2,2,4);imshow(Ib);
    
    end
    
    end
    
    function [ B ] = myrotate( A,O,angle )
    %定义我的旋转函数,表示A当中的点,绕O点,逆时针旋转angle角度后的点
    x0 = O(1);y0 = O(2);
    x1s = A(:,1);
    y1s = A(:,2);
    
     x1new = (x1s - x0)*cos(angle) - (y1s - y0)*sin(angle) + x0 ;
     y1new = (x1s - x0)*sin(angle) + (y1s - y0)*cos(angle) + y0 ;
     B = [x1new y1new];
    
    % myrotate([1 1],[2 1],-pi/4)
    end
    
    function [ Coordinates_normalized ] =  normalize( Coordinates )
    %此函数对于输入的二维坐标,按照x轴向映射的【0,1】区间
    %   对输入的坐标点归一化
    x = Coordinates(:,1);
    %y =  Coordinates(:,2);
    xmax = max(x);
    %xmin = min(x);%这个应该是等于0的
    Coordinates_normalized = Coordinates./xmax;
    end
    
    function pic_extract(video_name,freq)
    %视频中提取图像函数,video_name表示视频名称,freq表示每几帧提取一次
    
    %video_name = '1';
    VideoAd = VideoReader([video_name '.mp4']);%输入视频位置
    numFrames = VideoAd.NumberOfFrames;% 帧的总数
    videoF = VideoAd.FrameRate;%视频采集速率,每秒走几帧
    videoD = VideoAd.Duration; %时间
    %numname=3;%图片名称长度
    %nz = strcat('%0',num2str(numname),'d');
    T = 1*freq;%提取帧数间隔,这里设定每1秒提取一一帧
    %T=1;%提取帧数间隔,这里设定每1秒提取一一帧
    
    
    if ~exist(video_name,'dir')==1
       mkdir(video_name);%生成文件夹存放图片
    end
    i=1;
    for k = 1 :T: numFrames%? ? ?
        numframe = read(VideoAd,k);%读取第几帧
        %num=sprintf(nz,i);%i为保存图片的序号
        dir_path = strcat(video_name,'\',num2str(i));
        if ~exist(dir_path,'dir')==1
            mkdir(dir_path);%生成文件夹存放图片
        end
        imwrite(numframe,strcat(video_name,'\',num2str(i),'\',num2str(i),'.png'),'png'); %写入图片
        i=i+1;
    end
    
    end
    
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    %% 处理左图的程序
    clc
    clear
    addpath('../');%添加上层目录,方便调用函数
    D = dir('./*.bmp');
    D = struct2cell(D);%结构体不好调用,转成元胞数组来使用
    pics_name = D(1,:);
    for pic_name = pics_name %cpicName = 'a001.bmp'
        picName = pic_name{1};
        picN = strsplit(picName,'.');
        picN = picN{1};
        A = imread(picName);
        imshow(A);
        I = im2bw(A,100/255); %对图像二值化处理
        I = denoise(I);
        imshow(I);
        [ coord,L ] = edge_extract( ~I );
        imshow(L);
        x = coord(:,2);
        y = coord(:,1);
        inds = find(y == max(y));
        ind_max = max(inds);
        ind_min = min(inds);
        X0 = [max(y) (x(ind_min)+x(ind_max))/2 ];%使用矩阵坐标,即ij坐标
        point_num = size(coord,1);
        Coordinates = (repmat(X0,point_num,1) - coord);
        Coord =  normalize( Coordinates );%一个归一化处理
    
        imwrite(L,[picN 'E.png']); %写入图片
        csvwrite([picN '.csv'],Coord);
    
        figure;
        pic_std = plot(Coord(1:end,1),Coord(1:end,2),'r.');
        saveas(pic_std,[picN 'STD.jpg']); %写入图片
    end
    
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    %% 处理右图的程序
    clc
    clear
    addpath('../');%添加上层目录,方便调用函数
    D = dir('./*.bmp');
    D = struct2cell(D);%结构体不好调用,转成元胞数组来使用
    pics_name = D(1,:);
    for pic_name = pics_name %picName = 'b001.bmp'
        picName = pic_name{1};
        picN = strsplit(picName,'.');
        picN = picN{1};
        A = imread(picName);
        imshow(A);
        I = im2bw(A,60/255); %对图像二值化处理
        I = denoise(I);
        imshow(I);
        [ coord,L ] = edge_extract( ~I );%绕一圈儿,回到原点
        imshow(L);
        x = coord(:,2);
        y = coord(:,1);
        inds = find(y == max(y));
        ind_max = max(inds);
        ind_min = min(inds);
        X0 = [max(y) (x(ind_min)+x(ind_max))/2 ];%使用矩阵坐标,即ij坐标
    
        %找到上顶点的索引,做旋转矫正
        inds_up = find(y == min(y));
        ind_up_max = max(inds_up);
        ind_up_min = min(inds_up);
        Y0 = [min(y) (x(ind_up_min)+x(ind_up_max))/2 ];
        ang = -atan((X0(2)-Y0(2))/(X0(1)-Y0(1)));
        if abs(ang) > 0.01%角度太小就不用转了吧
        coord  = myrotate( coord,X0,ang);
        end
    
        %可视化矫正结果,调试用,只能调试到这一步时用,再往下X0不动的
    %     if min(min(coord))<1
    %         coord = coord-min(min(coord))+1;
    %     end
    %     [low_num,col_num] = size(I);
    %     L_test = sparse(round(coord(:,1)),round(coord(:,2)),1,low_num+100,col_num+100);
    %     L_test = full(L_test);
    %     imshow(L_test)
    %     figure(99);
    %     scatter(coord(:,1),coord(:,2))
    
        x = coord(:,2);
        y = coord(:,1);
    %     inds = find(y == max(y));
    %     ind_max = max(inds);
    %     ind_min = min(inds);
    %     X0 = [max(y) (x(ind_min)+x(ind_max))/2 ];%使用新选的X0
    
        inds_half = ((x<(X0(2)-1))&(y<=(X0(1))));
        coord_half = coord(inds_half,:);
    
        point_num = size(coord_half,1);
        coord = coord_half;
        Coordinates = (repmat(X0,point_num,1) - coord);
        % test 可视化结果
        %edge_rebuid = sparse(round(Coordinates(:,1))+10,round(Coordinates(:,2))+10,1,500,200);
        %min(round(Coordinates(:,2))+1)
        % imshow(full(edge_rebuid));
        % figure;
        % scatter(points_coord(:,1),points_coord(:,2))
        Coord =  normalize( Coordinates );%一个归一化处理
        imwrite(L,[picN 'E.png']); %写入图片
        csvwrite([picN '.csv'],Coord);
    
        figure;
        pic_std = plot(Coord(1:end,1),Coord(1:end,2),'r.');
        saveas(pic_std,[picN 'STD.jpg']); %写入图片
    
    %     pause(2);
    %     scatter(Coord(:,1),Coord(:,2));
    
    end
    
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    %% 主成分分析
    function [ coef ] = find_symbol_coef( exp,item )
    %输入一个表达式,即一个项,提取项前面的系数
    %   author:lusong
    %exp = Phi 很鲁棒
    charexp=char(exp);
    item=char(item);
    indexItem=strfind(charexp,item); %获取存在item项的指标
    indexPseudo=union(strfind(charexp,['(',item]),strfind(charexp,[item,'^'])); %获取伪指标
    indexItem=setdiff(indexItem,indexPseudo); %获取真正的item指标
    %循环计算各个item位置的系数
    itemLen=length(item);
    expLen=length(charexp);
    coef=sym(0);
    for i=1:length(indexItem)
        index=indexItem(i); %计算当前item项的位置
        cache=sym(1); %存储当前项的系数
        if index~=1 && charexp(index-1)=='*'
            indexFront=index-2; %初始化系数项的前指标
            while indexFront~=1 && charexp(indexFront-1)~=' '
                indexFront=indexFront-1;
            end
            cache=cache*sym(charexp(indexFront:index-2));
        end
        if index+itemLen~=expLen && charexp(index+1)=='*'
            indexBack=index+2; %初始化系数项的后指标
            while indexFront+itemLen~=expLen && charexp(indexBack+1)~=' '
                indexBack=indexBack+1;
            end
            cache=cache*sym(charexp(index+2:indexBack));
        end
        coef=coef+cache;
    end
    end
    
    clc;
    clear;
    close all;
    syms x y t;
    eta = 104;gamma = 20.9;rho = 964;
    theta_e = degtorad(53);alpha = degtorad(15);g = 9.8;
    num = 6;
    a = sym('a',[num,1]);
    adot = sym('adot',[num,1]);
    
    %% Phi的计算
    H = (x-a(1)).*(a(2)-x).*(a(3)+a(4)*x);
    Y = (x-a(1))^0.*(a(2)-x).^0*(a(5)+a(6)*x);
    h = H.*(1-(y./Y)^0.5);
    HY = H.*Y;
    HY_dot = conj(gradient(HY,a)')*adot;
    HY_dot = simplify(HY_dot);
    latex(HY_dot)
    V = -1/(HY).*int(HY_dot,x,a(1),x);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%这个积分算不出来,因为0.5次方的存在。
    Y_dot = conj(gradient(Y,a)')*adot;
    W = 1./Y*(Y_dot+V.*diff(Y,x));
    
    %Phi_intvar = 3*eta*(V.^2+(y*W).^2)./h;
    Phi_intvar = W.^2;
    %Phi_intvar = simplify(Phi_intvar);
    %latex(Phi_intvar)
    %Phi_intvar_resy = int(Phi_intvar,y,-Y,Y);
    %Phi = 0.5*(int(Phi_intvar_resy,x,a(1),a(2)));
    Phi_intvar_resx= int(Phi_intvar,x,a(1),a(2));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%这个积分算不出来,被积函数算不出来
    Phi = 0.5*(int(Phi_intvar_resx,y,-Y,Y));
    
    %% 提取Phi关于adot的系数
    A = a*conj(a');
    Phi = sum(sum(A.*(adot*conj(adot'))));
    
    adot_multi_matrix = adot*conj(adot');
    for i=1:num
        for j=1:num
            item = adot_multi_matrix(i,j);
            xi(i,j) = find_symbol_coef(Phi,item);
        end
    end
    xi = 2.*xi;
    
    % item = adot(2)*adot(1);
    %  find_symbol_coef(Phi,item)
    
    %% 势能函数
    int_var = 0.5.*gamma.*theta_e.^2+0.5.*gamma.*(diff(h,x).^2+diff(h,y).^2);
    remainder_term = 0.5.*rho.*g.*h^2.*sin(alpha)-rho.*g.*x.*h.*cos(alpha);
    A_temp = int(int_var,y,-Y,Y);
    A = int(A_temp,x,a(1),a(2))+remainder_term;%%%%%%%%%%积分也不好积
    
    %% 求解DOEs
    xi_inv = inv(xi);
    odefun = @(t,a) xi\gradient(A);%%%改
    a0 = ones(1,num);
    [T,a_cal] = ode45(odefun,[1,1],a0);
    
    %% 计算出来的结果以及画图
    H_cal = subs(H,[a1 a2 a3 a4],[a_cal(1) a_cal(2) a_cal(3) a_cal(4)]);
    Y_cal = subs(Y,[a1 a2 a5 a6],[a_cal(1) a_cal(2) a_cal(5) a_cal(6)]);
    H_cal_fun = matlabFunction(H_cal,x);
    xx = a(1):0.01:a(2);
    plot(xx,H_cal_fun(xx));
    
    
    clc
    clear
    close all
    format short
    digits(5)
    %% 设定维数
    p = 2;
    d = 2;
    %% 读入数据,保存
    D = dir('../data/*.csv');%%%%%%%%%%%%%%%%%%%%%%%%
    D = struct2cell(D);%结构体不好调用,转成元胞数组来使用
    csvs_name = D(1,:);
    N = 0;%这个表示数据的个数
    %csvs_name = csvs_name(5:6);
    for csv_name = csvs_name %cpicName = 'a001.bmp'
        N = N+1;
        csvName = csv_name{1};
        Data{N} = load(['../data/' csvName]);%%%%%%%%%%%%%%%%%%%%%
    end
    %clear csv_name csvN csvName csvs_name D
    %% 定义原空间基,生成模型函数,由模型函数通过拟合得到点的坐标表示
    %p = 15;%原空间维数
    %%根据p自动生成模型函数
    eval_poly = 'a(1)';
    for i = 2:p
        eval_poly = [eval_poly '+a(' num2str(i) ').*x.^' num2str(i-1)];
    end
    eval_poly = ['modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(' eval_poly ');'];
    eval(eval_poly);
    
    %modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*sin(x)+a(3).*cos(x)+a(4).*1./x+a(5).*x+a(6).*x.^2);
    %   modelfun = @(a,x)(x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*x+a(3).*x.^2+...
    %       a(4).*x.^3))%+a(5).*x.^4+a(6)*x.^5))+...
    %       a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
    %         modelfun = @(a,x)((a(1)+a(2).*x.^2+a(3).*x.^2+...
    %         a(4).*x.^3+a(5).*x.^4+a(6)*x.^5+...
    %         a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
    
    %%做一波拟合,求出系数,进行拟合,拟合结果为A
    %%A中存的是数据点的坐标表示
    for k = 1:N
        data = Data{k};
        data(:,2) = abs(data(:,2));%取绝对值
        x = data(:,1);
        y = data(:,2);
        %y = y./(x.*(1-x));%%%%%%%%如果将x(1-x)项挪到左边,拟合会出问题
        a0 = ones(1,p);
        a = nlinfit(x,y,modelfun,a0);
        A(:,k) = a';
    end
    %clearvars -except A modelfun x y
    
    %% 开始奇异值分解降维,得到新空间的原点和坐标轴,坐标表示
    %%由奇异值分解做主成分分析,A的每一列为原空间下的系数,U的每一列为主成分
    [p,N] = size(A);
    %d = 5;%新空间维数
    x0 = mean(A,2);%计算样本均值
    A_shift = A - repmat(x0,1,N);%中心平移每个训练样本,A_shift为原数据点的平移表示
    
    %%%%对于A_shift的补充处理
    syms x;
    base_vec_T = x.^(0.5).*(1-x).^(0.5).*x.^[0:p-1];
    base_vec = conj((base_vec_T))';
    W = int(base_vec*base_vec_T,[0,1]);
    L_prime = chol(W);
    L_prime_A_shift = L_prime*A_shift;
    
    %%%%
    [U,S,V] = svd(L_prime_A_shift);
    B = L_prime\U;
    if size(S,2) == 1
        S_diag = S(1);
    else
        S_diag = diag(S);
    end
    R = cumsum(S_diag)./sum(S_diag);
    Coord_new = B(:,1:d);
    
    
    %% 新空间的坐标轴和原点的函数表示
    %%生成新的空间中的基函数,找到d个新空间基函数及原点坐标
    %%Coord_new x0中存的是新坐标系,fi存的是新的函数空间的原点和坐标轴
    for k = 1:d;%取出第k个基函数
        syms x;
        c = B(:,k);
        exp = modelfun(c,x);
        exp = simplify(exp);%第k个基函数的表达式
        f{k} = exp;
        digits(2);
        latex(vpa(exp,2));%写成latex表达式黏贴到解释器中
    end
    
    f0 = modelfun(x0,x);
    f0 = simplify(f0);
    latex(vpa(f0,2));
    % for i = 1:d
    %     f_fun{i} = matlabFunction(f{i});
    % end
    % f0_fun = matlabFunction(f0);
    
    %% 求原坐标系中的投影后在原空间中的位置
    
    
    
    %A_new_coord = Coord_new'*A_shift;%A_new_coord表示原数据点在新空间中的表示
    
    syms x;
    for k = 1:N
        a = A_shift(:,k);
        pk_f = modelfun(a,x);
        for j = 1:d
            A_new_coord(j,k) = int(pk_f.*f{j},0,1);
        end
    end
    
    %%函数运算
    for k = 1:N
        A_shadow_f = f0;
        for i = 1:d
            A_shadow_f = A_shadow_f+f{i}*A_new_coord(i,k);
        end
        A_shadow_f = simplify(A_shadow_f);
        A_shadow_fs{k} = A_shadow_f;
        A_shadow_fs_fun{k} = matlabFunction(A_shadow_f);
    end
    
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % %%坐标运算,可以注释掉
    % A_shadow_coord = Coord_new*A_new_coord+repmat(x0,1,N);%A_shadow_coord表示投影坐标点在原坐标系下的表示
    % syms x;
    % for k = 1:N
    %     a = A_shadow_coord(:,k);
    %     A_shadow_coord_f = modelfun(a,x);
    %     A_shadow_coord_fs{k} = simplify(A_shadow_coord_f);
    % end
    %
    % %%不考虑舍入误差的情况下,可以证明坐标运算和函数运算的得到的函数表达式是一样的,验证代码如下
    % k=5;
    % diff = A_shadow_coord_fs{k}-A_shadow_fs{k}
    % diff = simplify(diff)
    % digits(2);
    % latex(vpa(diff,2))
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %% 画图和求误差
    Errors = [];
    NN = 1000;
    h = 1/NN;
    xx = linspace(0,1,NN+1);
    close all;
    for k = 1:N
        %%原数据点
        data = Data{k};
        data(:,2) = abs(data(:,2));%取绝对值
        x = data(:,1);
        y = data(:,2);
        %%拟合曲线
        a = A(:,k);
        data_fit_f = @(x) modelfun(a,x);
        %%投影函数
        A_shadow_f_fun = A_shadow_fs_fun{k};
        %%原点函数
        f0_fun = matlabFunction(f0);
    
        %%开始画图
        figure(k);
        scatter(x,y,'.','MarkerEdgeColor','b',...
                  'MarkerFaceColor','c',...
                  'LineWidth',0.5);
        hold on;
    
        y_fit = data_fit_f(xx);
        y_shadow = A_shadow_f_fun(xx);
        y_f0 = f0_fun(xx);
        plot(xx,y_fit,'g');
        plot(xx,y_shadow,'black');
        plot(xx,y_f0,'r');
        legend('数据','拟合','投影','原点');
    
        %%求拟合曲线和投影函数之间的L2误差
        %Errors(end+1) = norm((y_fit-y_shadow),2);
        Errors(end+1) = sum(abs(y_fit-y_shadow)*h);
    end
    
    Errors_mean = mean(Errors)
    
    
    
    %% test
    % close all;
    % plot(xx,xx);
    % hold on;
    % plot(xx,xx.^9)
    
    clc
    clear
    close all
    format short
    digits(5)
    %% 设定维数
    p = 2;
    d = 1;
    %% 读入数据,保存
    D = dir('../data2/*.csv');%%%%%%%%%%%%%%%%%%%%%%%%
    D = struct2cell(D);%结构体不好调用,转成元胞数组来使用
    csvs_name = D(1,:);
    N = 0;%这个表示数据的个数
    %csvs_name = csvs_name(5:6);
    for csv_name = csvs_name %cpicName = 'a001.bmp'
        N = N+1;
        csvName = csv_name{1};
        Data{N} = load(['../data2/' csvName]);%%%%%%%%%%%%%%%%%%%%%
    end
    %clear csv_name csvN csvName csvs_name D
    %% 定义原空间基,生成模型函数,由模型函数通过拟合得到点的坐标表示
    %p = 15;%原空间维数
    %%根据p自动生成模型函数
    eval_poly = 'a(1)';
    for i = 2:p
        eval_poly = [eval_poly '+a(' num2str(i) ').*x.^' num2str(i-1)];
    end
    eval_poly = ['modelfun = @(a,x) x.^(1).*(1-x).^(1).*(' eval_poly ');'];
    eval(eval_poly);
    
    %modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*sin(x)+a(3).*cos(x)+a(4).*1./x+a(5).*x+a(6).*x.^2);
    %   modelfun = @(a,x)(x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*x+a(3).*x.^2+...
    %       a(4).*x.^3))%+a(5).*x.^4+a(6)*x.^5))+...
    %       a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
    %         modelfun = @(a,x)((a(1)+a(2).*x.^2+a(3).*x.^2+...
    %         a(4).*x.^3+a(5).*x.^4+a(6)*x.^5+...
    %         a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
    
    %%做一波拟合,求出系数,进行拟合,拟合结果为A
    %%A中存的是数据点的坐标表示
    for k = 1:N
        data = Data{k};
        data(:,2) = abs(data(:,2));%取绝对值
        x = data(:,1);
        y = data(:,2);
        %y = y./(x.*(1-x));%%%%%%%%如果将x(1-x)项挪到左边,拟合会出问题
        a0 = ones(1,p);
        a = nlinfit(x,y,modelfun,a0);
        A(:,k) = a';
    end
    %clearvars -except A modelfun x y
    
    %% 开始奇异值分解降维,得到新空间的原点和坐标轴,坐标表示
    %%由奇异值分解做主成分分析,A的每一列为原空间下的系数,U的每一列为主成分
    [p,N] = size(A);
    %d = 5;%新空间维数
    x0 = mean(A,2);%计算样本均值
    A_shift = A - repmat(x0,1,N);%中心平移每个训练样本,A_shift为原数据点的平移表示
    
    %%%%对于A_shift的补充处理
    syms x;
    base_vec_T = x.^(1).*(1-x).^(1).*x.^[0:p-1];
    base_vec = conj((base_vec_T))';
    W = int(base_vec*base_vec_T,[0,1]);
    L_prime = chol(W);
    L_prime_A_shift = L_prime*A_shift;
    
    %%%%
    [U,S,V] = svd(L_prime_A_shift);
    B = L_prime\U;
    if size(S,2) == 1
        S_diag = S(1);
    else
        S_diag = diag(S);
    end
    R = cumsum(S_diag)./sum(S_diag);
    Coord_new = B(:,1:d);
    
    
    %% 新空间的坐标轴和原点的函数表示
    %%生成新的空间中的基函数,找到d个新空间基函数及原点坐标
    %%Coord_new x0中存的是新坐标系,fi存的是新的函数空间的原点和坐标轴
    for k = 1:d;%取出第k个基函数
        syms x;
        c = B(:,k);
        exp = modelfun(c,x);
        exp = simplify(exp);%第k个基函数的表达式
        f{k} = exp;
        digits(2);
        latex(vpa(exp,2));%写成latex表达式黏贴到解释器中
    end
    
    f0 = modelfun(x0,x);
    f0 = simplify(f0);
    latex(vpa(f0,2));
    % for i = 1:d
    %     f_fun{i} = matlabFunction(f{i});
    % end
    % f0_fun = matlabFunction(f0);
    
    %% 求原坐标系中的投影后在原空间中的位置
    
    
    
    %A_new_coord = Coord_new'*A_shift;%A_new_coord表示原数据点在新空间中的表示
    
    syms x;
    for k = 1:N
        a = A_shift(:,k);
        pk_f = modelfun(a,x);
        for j = 1:d
            A_new_coord(j,k) = int(pk_f.*f{j},0,1);
        end
    end
    
    %%函数运算
    for k = 1:N
        A_shadow_f = f0;
        for i = 1:d
            A_shadow_f = A_shadow_f+f{i}*A_new_coord(i,k);
        end
        A_shadow_f = simplify(A_shadow_f);
        A_shadow_fs{k} = A_shadow_f;
        A_shadow_fs_fun{k} = matlabFunction(A_shadow_f);
    end
    
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % %%坐标运算,可以注释掉
    % A_shadow_coord = Coord_new*A_new_coord+repmat(x0,1,N);%A_shadow_coord表示投影坐标点在原坐标系下的表示
    % syms x;
    % for k = 1:N
    %     a = A_shadow_coord(:,k);
    %     A_shadow_coord_f = modelfun(a,x);
    %     A_shadow_coord_fs{k} = simplify(A_shadow_coord_f);
    % end
    %
    % %%不考虑舍入误差的情况下,可以证明坐标运算和函数运算的得到的函数表达式是一样的,验证代码如下
    % k=5;
    % diff = A_shadow_coord_fs{k}-A_shadow_fs{k}
    % diff = simplify(diff)
    % digits(2);
    % latex(vpa(diff,2))
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %% 画图和求误差
    Errors = [];
    NN = 1000;
    h = 1/NN;
    xx = linspace(0,1,NN+1);
    close all;
    for k = 1:N
        %%原数据点
        data = Data{k};
        data(:,2) = abs(data(:,2));%取绝对值
        x = data(:,1);
        y = data(:,2);
        %%拟合曲线
        a = A(:,k);
        data_fit_f = @(x) modelfun(a,x);
        %%投影函数
        A_shadow_f_fun = A_shadow_fs_fun{k};
        %%原点函数
        f0_fun = matlabFunction(f0);
    
        %%开始画图
        figure(k);
        scatter(x,y,'.','MarkerEdgeColor','b',...
                  'MarkerFaceColor','c',...
                  'LineWidth',0.5);
        hold on;
    
        y_fit = data_fit_f(xx);
        y_shadow = A_shadow_f_fun(xx);
        y_f0 = f0_fun(xx);
        plot(xx,y_fit,'g');
        plot(xx,y_shadow,'black');
        plot(xx,y_f0,'r');
        legend('数据','拟合','投影','原点');
    
        %%求拟合曲线和投影函数之间的L2误差
        %Errors(end+1) = norm((y_fit-y_shadow),2);
        Errors(end+1) = sum(abs(y_fit-y_shadow)*h);
    end
    
    Errors_mean = mean(Errors)
    
    
    
    %% test
    % close all;
    % plot(xx,xx);
    % hold on;
    % plot(xx,xx.^9)
    
    \begin{lstlisting}[language=Matlab,frame=shadowbox]
    %% 最小二乘
    clc
    clear
    D = dir('data/*.csv');
    D = struct2cell(D);
    names = D(1,:);
    data = [];
    for name = names
    name = name{1};
    data = [data;load(['data/' name])];
    end
    x = data(:,1);
    y = data(:,2);
    %modelfun = @(a,x)(x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*x+a(6)*x*4+a(3).*x.^2+a(4).*x.^3+a(5).*x.^5))%+a(6)*x.*4));
    modelfun = @(a,x)(x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*x.*100+a(3).*x.^2+a(4).*x.^3+a(5).*x.^4+a(6)*x.^5));
    a0 = ones(1,6)*100;
    modelfun(a0,0)
    modelfun(a0,1)
    a = nlinfit(x,y,modelfun,a0)
    a(abs(a)<2) = 0
    plot(x,modelfun(a,x))
    
    
    展开全文
  • 二维数据主成分分析的原理—— PCA投影PCA与线性回归的区别matlab自带PCA函数代码 数据降维 二维数据->一维数据、 降维到一维之后的数据: 三维数据->二维数据 上图中的数据可以投影到二维平面上,大致...

  • 主成分分析的原理—— PCA投影
  • PCA与线性回归的区别
  • matlab自带PCA函数
  • 代码
  • 主成分回归

    判定自变量相关程度大小,即变量间是否存在严重的多重共线性。若存在就用主成分分析,将存在多重共线性的变量合并为一个新的变量,然后再和其它自变量纳入回归。这整个过程叫主成分回归。

    数据降维

    二维数据->一维数据、


    降维到一维之后的数据:


    三维数据->二维数据


    上图中的数据可以投影到二维平面上,大致在如下平面上:


    投影之后如下图所示:


    主成分分析的原理—— PCA投影

    假设有一组二维数据如下图所示:


    若将上图中的数据投影到图中的红色直线上,注意到原始数据到投影数据的距离为蓝色线段,且距离最短。PCA求解的目标是:将原始数据投影到低维平面上,使得图中蓝色线段的和达到最小,即投影误差最小。


    换一个投影方向,投影到图中玫红色的直线上,投影误差达到最大化。


    将二维数据投影成一维数据,需要寻找一个投影方向,或者是一个n维向量


    将三维数据降为二维数据,需要寻找两个n维向量,组成一个投影平面:



    PCA与线性回归的区别

    在PCA由二维投影成一维的过程中,与线性回归类似,都是寻找一条直线,但是本质却不同。左图是线性回归中的预测直线,预测值与真实值之间的误差是图中蓝色线段,可见蓝色线短和y轴平行;右图是PCA的投影方向,蓝色线段是投影误差,可见蓝色线段垂直于投影方向,是最短距离。除此之外,线性回归需要使用标记数据,获取样本的类别标记,而PCA算法则是无监督的,不需要使用样本的类别标记。


    matlab自带PCA函数

    数据集X(每行为一个样本,行数为样本数)

    • coeff = pca(X)
    • coeff = pca(X,Name,Value)
    • [coeff,score,latent] = pca(___)
    • [coeff,score,latent,tsquared] = pca(___)
    • [coeff,score,latent,tsquared,explained,mu] = pca(___)

    Input Argument 0
    X :--数据集 假设n个样本, 每个样本p维,则 X是n-by-p的matrix

    Input Argument 1
    'Algorithm' — Principal component algorithm
    'svd' (default) | 'eig' | 'als'
    解释:PCA 涉及到求协方差矩阵的特征向量, 在matlab 有三种算法
    默认 :SVD,
    eig (Eigenvalue decomposition )算法, 此算法当n(number of examples) > p (features) 时,速度快于SVD,但是计算的结果没有SVD精确
    als( Alternating least squares )算法,此算法为了处理数据集X中有少许缺失数据时的情况(i.e 0), 但是对于X为稀疏数据集(缺失数据过多)时,不好用

    Input Argument 2
    'Centered' — Indicator for centering columns
    true (default) | false
    解释:选择是否对数据进行中心化, 也是数据的特征是否进行零均值化(i.e.按列减去均值, 为了得到covariance matrix), 如果选择了true,  则可用score*coeff'恢复中心化后的X, 若选择了false,则可用score*coeff'恢复原始的X

    Input Argument 3
    'Economy' — Indicator for economy size output
    true (default) | false
    解释: 有时候输出的coeff(映射矩阵p-by-p)过大, 而且是没有必要的(因为我们要降维),所以可以只输出coeff(以及score,latent)的前d列,
    d是数据集的自由度,数据集没NAN的时候d=n-1; 具体的解释见matlab.总之如果将看见完整的PCA结果,可以设置为false.
    默认:true ,(默认ture以后对于初次使用matlab这个函数的人非常迷惑).

    Input Argument 4
    'NumComponents' — Number of components requested
    number of variables (default) | scalar integer
    解释:输出指定的components 也就是更为灵活的Economy,但是经过试验发现指定成分数 仅在小于d(自由度)时有效,大于d时无效;
    默认: number of variables ( i.e p,特征数目)

    Input Argument 5
    'Rows' — Action to take for NaN values
    'complete' (default) | 'pairwise' | 'all'
    解释: 此选项是为了智能处理数据集X中含有NAN的情况
    complete: 计算之前.移除X中含有NAN的行(i.e 样本),计算完成后,含有NAN的行被重新插入到score and tsquared相应的位置中(coeff呢?)
    pairwise : 首先这个选项必须配合 'Argorithm'中 'eig'进行使用.如果没有指定'eig'(默认svd),当指定pairwise时,则会自动切换为eig; 指定为svd,则会发送warning message,然后自动切换为eig;若指定为als, 则会发送warning message然后忽略 'Rows'此选项。 成功使用此选项时,若计算协方差(i,j)处值时遇到NAN,则使用X中第i或j列中不含NAN的行此处来计算的协方差值.
    all : 当确定X中无缺失数据,使用'all',则pca会使用X中所有的数据,当遇到NAN时则会自动终止.

    Input Argument 6
    'Weights' — Observation weights
    ones (default) | row vector
    解释: 基于observations(i.e 样本)的权重pca,有需求的可以自己查查

    Input Argument 7
    'VariableWeights' — Variable weights
    row vector | 'variance'
    解释:基于variables(i.e.features)的权重pca,有需求的自己查
    默认: 无默认值, 也就是默认不使用此选项

    Input Argument 8
    'Coeff0' — Initial value for coefficients
    matrix of random values (default) | p-by-k matrix
    解释: Initial value for the coefficient matrix coeff, 不是太看的懂,但是要配合'Algorithm'中的'als'使用
    默认:   p-by-random

    Input Argument 9
    'Score0' — Initial value for scores
    matrix of random values (default) | k-by-m matrix
    解释: Initial value for scores matri score.不是太看的懂,但是要配合'Algorithm'中的 'als'使用
    默认 : n-by-random

    Input Argument 10
    'Options' — Options for iterations
    structure(此用法是个结构体)
    解释:用于迭代的选项,仅配合'Algorithm'中的'als'使用. 因为'als'是使用迭代的方法进行计算的、对这个不感兴趣, 有兴趣的可以去help一下附上help中的使用方法 opt = statset('pca'); opt.MaxIter = 2000; coeff =pca(X,'Options',opt);

    Output Argument 1
    coeff : 主成分系数 应该就是协方差矩阵的特征向量矩阵(也就是映射矩阵).完整输出的情况下是一个p-by-p 的matrix.每一列都是一个特征向量.按对应的特征值的大小,从大到小进行排列.

    Output Argument 2
    score: 进行旋转(也就是利用映射矩阵coeff进行)后的结果i.e. score = X * coeff. n-by-p matrix这里有个坑 如果你使用pca时使用的是默认的中心化(i.e 不对’Centered’设置’false’),拿X * coeff 和score对比的时候,记得把X中心化后再乘以coeff,之后再和score对比…;同样如果pca使用的是默认值, 恢复的X = score * coeff’ (注意转置)是中心化后的数据

    Output Argument 3
    latent: 主成分方差 也就是各特征向量对应的特征值,从大到小进行排列

    Output Argument 4
    tsquared :层次不够 无法解释......

    Output Argument 5
    explained : 每一个主成分所贡献的比例,可以更直观的选择所需要降维的维数了,不用再用特征值去求了

    Output Argument 6
    mu: X 按列的均值,当前仅当 'Centered'置于'true'(默认值)时才会返回此变量
     

    代码

    clc;
    clear;
    
    load fisheriris
    data=meas;
    
    % data行为不同样本的特征,列为不同特征,latent长度与data列数相同,res为维度累计占比
    % coeff按对应的特征值的大小,从大到小进行排列.
    [coeff,score,latent,tsquare] = pca(data);%我们这里需要他的pc和latent值做分析
    % 根据res结果决定降到final维
    res=cumsum(latent)./sum(latent);
    final=2;
    % 取coeff中的1:final列来做最后的变换矩阵
    data_after_PCA=score(:,1:final);
    % 散点图表示每个样例在pc1和pc2两个成分上的位置,可以看出大量的鸢尾花数据基本分为两类,每个点表示一个鸢尾花数据
    scatter(data_after_PCA(:,1),data_after_PCA(:,2))
    
    展开全文
  • 上一篇中我们详细介绍推导了主成分分析法的原理,并基于Python通过自编函数实现了挑选主成分的过程,而在Python与R中都有比较成熟的主成分分析函数,本篇我们就对这些方法进行介绍: R 在R的基础函数中就有...

    上一篇中我们详细介绍推导了主成分分析法的原理,并基于Python通过自编函数实现了挑选主成分的过程,而在Python与R中都有比较成熟的主成分分析函数,本篇我们就对这些方法进行介绍:

     

    R

    在R的基础函数中就有主成分分析法的实现函数princomp(),其主要参数如下:

    data:要进行主成分分析的目标数据集,数据框形式,行代表样本,列代表变量

    cor:逻辑型变量,控制是否使用相关系数进行主成分分析

    scores:逻辑型变量,控制是否计算每个主成分的得分

    我们使用了R中自带的数据集USJudgeRating来进行演示,这是一个包含43个样本,12个连续型实自变量的数据集,适合来演示PCA,这里我们在其自带方法的基础上,使用自编函数来对训练后的数据进行一步到位的,基于用户自定义贡献率阈值(默认0.8)来生成所需主成分值,具体过程如下:

    > rm(list=ls())
    > 
    > #载入律师评价数据
    > data(USJudgeRatings)
    > data <- USJudgeRatings
    > 
    > #对律师评价数据进行主成分分析,这里设置使用相关系数进行主成分分析
    > data.pr <- princomp(data,cor=T,scores=T)
    > 
    > #自定义目标主成分值生成函数
    > My_Choose <- function(data.pr,a=0.8){
    +   contribute <- data.pr$sdev^2/sum(data.pr$sdev^2)
    +   c <- 0
    +   i <- 0
    +   #根据累计贡献率挑选所需最少的主成分
    +   while(c<=a){
    +     i <- i + 1
    +     c <- c + contribute[i]
    +   }
    +   cat('达到要求,已采纳的主成分累计贡献率达到:','\n',round(c,4),'\n')
    +   return(predict(data.pr)[,1:i])
    + }
    > 
    > #显示主成分分析的结果,以及因子载荷
    > summary(data.pr, loadings = T)
    Importance of components:
                              Comp.1     Comp.2    Comp.3
    Standard deviation     3.1833165 1.05078398 0.5769763
    Proportion of Variance 0.8444586 0.09201225 0.0277418
    Cumulative Proportion  0.8444586 0.93647089 0.9642127
                               Comp.4      Comp.5      Comp.6
    Standard deviation     0.50383231 0.290607615 0.193095982
    Proportion of Variance 0.02115392 0.007037732 0.003107172
    Cumulative Proportion  0.98536661 0.992404341 0.995511513
                                Comp.7      Comp.8       Comp.9
    Standard deviation     0.140295449 0.124158319 0.0885069038
    Proportion of Variance 0.001640234 0.001284607 0.0006527893
    Cumulative Proportion  0.997151747 0.998436354 0.9990891437
                                Comp.10      Comp.11      Comp.12
    Standard deviation     0.0749114592 0.0570804224 0.0453913429
    Proportion of Variance 0.0004676439 0.0002715146 0.0001716978
    Cumulative Proportion  0.9995567876 0.9998283022 1.0000000000
    
    Loadings:
         Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
    CONT         0.933 -0.335                                   
    INTG -0.289 -0.182 -0.549 -0.174         0.370 -0.450  0.334
    DMNR -0.287 -0.198 -0.556  0.124  0.229 -0.395  0.467 -0.247
    DILG -0.304         0.164 -0.321  0.302  0.599  0.210 -0.355
    CFMG -0.303  0.168  0.207         0.448         0.247  0.714
    DECI -0.302  0.128  0.298         0.424 -0.393 -0.536 -0.302
    PREP -0.309         0.152 -0.214 -0.203         0.335 -0.154
    FAMI -0.307         0.195 -0.201 -0.507 -0.102              
    ORAL -0.313                      -0.246 -0.150              
    WRIT -0.311               -0.137 -0.306 -0.238         0.126
    PHYS -0.281         0.154  0.841 -0.118  0.299              
    RTEN -0.310        -0.173  0.184               -0.256 -0.221
         Comp.9 Comp.10 Comp.11 Comp.12
    CONT                               
    INTG -0.275 -0.109  -0.113         
    DMNR -0.199          0.134         
    DILG         0.383                 
    CFMG  0.143          0.166         
    DECI -0.258         -0.128         
    PREP -0.109 -0.680  -0.319   0.273 
    FAMI -0.223          0.573  -0.422 
    ORAL  0.300  0.256  -0.639  -0.494 
    WRIT         0.475           0.696 
    PHYS -0.266                        
    RTEN  0.756 -0.250   0.286         
    > 
    > #保存主成分数据
    > (pca_data <- My_Choose(data.pr,a=0.9))
    达到要求,已采纳的主成分累计贡献率达到: 
     0.9365 
                         Comp.1       Comp.2
    AARONSON,L.H.    0.59268305 -1.834557041
    ALEXANDER,J.M.  -2.40817832 -0.878164524
    ARMENTANO,A.J.  -0.22780349 -0.302316250
    BERDON,R.I.     -3.66071616 -0.596342199
    BRACKEN,J.J.     6.95244952  0.137356050
    BURNS,E.B.      -2.47442074 -1.422422886
    CALLAHAN,R.J.   -3.93743704  3.147711145
    COHEN,S.S.       8.09307094 -0.271058974
    DALY,J.J.       -3.70295623 -0.146659161
    DANNEHY,J.F.    -1.06028925  1.149046170
    DEAN,H.H.        0.34108279 -0.470502584
    DEVITA,H.J.      1.45772527 -1.232794755
    DRISCOLL,P.J.    0.67688211 -1.297855808
    GRILLO,A.E.      3.26409534 -0.566023927
    HADDEN,W.L.JR.  -1.36896588 -0.848756359
    HAMILL,E.C.      0.45009339 -0.182461748
    HEALEY.A.H.      2.93209773  0.383887639
    HULL,T.C.        0.72643742  0.463364335
    LEVINE,I.       -0.64887947  0.904423859
    LEVISTER,R.L.    4.40757580  2.444558329
    MARTIN,L.F.      1.84009523 -0.801922043
    MCGRATH,J.F.     3.15606059  0.118968056
    MIGNONE,A.F.     6.42980257 -1.452770363
    MISSAL,H.M.      0.08315548 -1.483307990
    MULVEY,H.M.     -3.39121775  0.132351096
    NARUK,H.J.      -4.60271493  0.479764153
    O'BRIEN,F.J.    -1.51784742 -0.427971456
    O'SULLIVAN,T.J. -3.51102566  0.001182987
    PASKEY,L.       -1.91559427  0.246170053
    RUBINOW,J.E.    -4.83537434 -0.448069023
    SADEN.G.A.      -1.11213340 -0.344165889
    SATANIELLO,A.G. -0.54875827  1.022509509
    SHEA,D.M.       -2.61362018 -0.399509080
    SHEA,J.F.JR.    -3.63917622 -0.162216363
    SIDOR,W.J.       6.94389214  0.330376531
    SPEZIALE,J.A.   -2.04961153  1.193213351
    SPONZO,M.J.     -1.22081130 -0.531924504
    STAPLETON,J.F.  -0.72842414 -0.914724377
    TESTO,R.J.       2.13120556  0.949624029
    TIERNEY,W.L.JR. -1.36785124  1.100438878
    WALL,R.A.        2.61916318  1.926895688
    WRIGHT,D.B.     -1.48026785 -0.556116054
    ZARRILLI,K.J.    0.92650698  1.440771500

    在得到累计贡献率高达0.9365的两个主成分之后,我们将主成分降维前后的数据的相关系数矩阵进行比较:

    > #降维前
    > cor(data)
                CONT       INTG       DMNR      DILG      CFMG
    CONT  1.00000000 -0.1331909 -0.1536885 0.0123920 0.1369123
    INTG -0.13319089  1.0000000  0.9646153 0.8715111 0.8140858
    DMNR -0.15368853  0.9646153  1.0000000 0.8368510 0.8133582
    DILG  0.01239200  0.8715111  0.8368510 1.0000000 0.9587988
    CFMG  0.13691230  0.8140858  0.8133582 0.9587988 1.0000000
    DECI  0.08653823  0.8028464  0.8041168 0.9561661 0.9811359
    PREP  0.01146921  0.8777965  0.8558175 0.9785684 0.9579140
    FAMI -0.02563656  0.8688580  0.8412415 0.9573634 0.9354684
    ORAL -0.01199681  0.9113992  0.9067729 0.9544758 0.9505657
    WRIT -0.04381025  0.9088347  0.8930611 0.9592503 0.9422470
    PHYS  0.05424827  0.7419360  0.7886804 0.8129211 0.8794874
    RTEN -0.03364343  0.9372632  0.9437002 0.9299652 0.9270827
               DECI       PREP        FAMI        ORAL
    CONT 0.08653823 0.01146921 -0.02563656 -0.01199681
    INTG 0.80284636 0.87779650  0.86885798  0.91139915
    DMNR 0.80411683 0.85581749  0.84124150  0.90677295
    DILG 0.95616608 0.97856839  0.95736345  0.95447583
    CFMG 0.98113590 0.95791402  0.93546838  0.95056567
    DECI 1.00000000 0.95708831  0.94280452  0.94825640
    PREP 0.95708831 1.00000000  0.98986345  0.98310045
    FAMI 0.94280452 0.98986345  1.00000000  0.98133905
    ORAL 0.94825640 0.98310045  0.98133905  1.00000000
    WRIT 0.94610093 0.98679918  0.99069557  0.99342943
    PHYS 0.87176277 0.84867350  0.84374436  0.89116392
    RTEN 0.92499241 0.95029259  0.94164495  0.98213227
                WRIT       PHYS        RTEN
    CONT -0.04381025 0.05424827 -0.03364343
    INTG  0.90883469 0.74193597  0.93726315
    DMNR  0.89306109 0.78868038  0.94370017
    DILG  0.95925032 0.81292115  0.92996523
    CFMG  0.94224697 0.87948744  0.92708271
    DECI  0.94610093 0.87176277  0.92499241
    PREP  0.98679918 0.84867350  0.95029259
    FAMI  0.99069557 0.84374436  0.94164495
    ORAL  0.99342943 0.89116392  0.98213227
    WRIT  1.00000000 0.85594002  0.96755639
    PHYS  0.85594002 1.00000000  0.90654782
    RTEN  0.96755639 0.90654782  1.00000000
    > #降维后
    > cor(pca_data)
                 Comp.1       Comp.2
    Comp.1 1.000000e+00 2.161808e-16
    Comp.2 2.161808e-16 1.000000e+00

    可以看出,降维后两个主成分之间的相关系数非常低,可以说它们几乎正交,说明主成分的结果非常有效:

     

    Python

    我们使用sklearn.decomposition中的PCA来实现主成分降维,其主要参数如下:

    n_components:这个参数可以帮我们指定希望PCA降维后的特征维度数目。最常用的做法是直接指定降维到的维度数目,此时n_components是一个大于等于1的整数。当然,我们也可以指定主成分的累计贡献率阈值,让PCA类自己去根据样本特征方差来决定降维到的维度数,此时n_components是一个(0,1]之间的数。当然,我们还可以将参数设置为"mle", 此时PCA类会用MLE算法根据特征的方差分布情况自己去选择一定数量的主成分特征来降维。我们也可以用默认值,即不输入n_components,此时n_components=min(样本数,特征数)

    whiten :判断是否进行白化。所谓白化,就是对降维后的数据的每个特征进行归一化,让方差都为1。对于PCA降维本身来说,一般不需要白化。如果你PCA降维后有后续的数据处理,可以考虑白化。默认值是False,即不进行白化。

    svd_solver:即指定奇异值分解SVD的方法,由于特征分解是奇异值分解SVD的一个特例,一般的PCA库都是基于SVD实现的。有4个可以选择的值:{‘auto’, ‘full’, ‘arpack’, ‘randomized’}。randomized一般适用于数据量大,数据维度多同时主成分数目比例又较低的PCA降维,它使用了一些加快SVD的随机算法。 full则是传统意义上的SVD,使用了scipy库对应的实现。arpack和randomized的适用场景类似,区别是randomized使用的是scikit-learn自己的SVD实现,而arpack直接使用了scipy库的sparse SVD实现。默认是auto,即PCA类会自己去在前面讲到的三种算法里面去权衡,选择一个合适的SVD算法来降维。一般来说,使用默认值就够了。

      除了这些输入参数外,有两个PCA类的输出值得关注。第一个是explained_variance_,它代表降维后的各主成分的方差值。方差值越大,则说明越是重要的主成分。第二个是explained_variance_ratio_,它代表降维后的各主成分的方差值占总方差值的比例,这个比例越大,则越是重要的主成分。

    我们选用datasets中自带的wine数据集作为演示数据,关于这个数据集可以参考前一篇的介绍,具体过程如下:

    from sklearn.decomposition import PCA
    from sklearn import datasets
    import numpy as np
    
    '''载入数据'''
    X,y = datasets.load_wine(return_X_y=True)
    
    '''显示待转换数据的形状'''
    print(X.shape)
    
    '''初始化PCA模型,这里选择希望从13个原始变量中产出三个主成分'''
    pca = PCA(n_components=3)
    
    '''将X导入设定好的模型中'''
    pca.fit(X)
    
    '''打印产出的主成分对应的方差贡献率'''
    print(pca.explained_variance_ratio_)
    
    '''打印产出的主成分对应的方差'''
    print(pca.explained_variance_)
    
    '''将模型产出的目标主成分保存'''
    X_new = pca.fit_transform(X)

    原始数据的形状以及得到的主成分的各自的方差贡献率、方差如下:

    下面计算原始数据的相关系数矩阵中元素的平均值与得到的主成分进行对比:

    '''计算原始数据相关系数矩阵的平均值'''
    print('原始相关系数矩阵元素的平均值:'+'\n'+str(np.mean(np.corrcoef(X))))
    
    '''计算保留的主成分之间的相关系数'''
    print('主成分的相关系数矩阵:'+'\n'+str(np.corrcoef(X_new.T)))
    
    '''计算主成分相关系数矩阵的平均值'''
    print('得到的主成分的相关系数矩阵元素的平均值:'+'\n'+str(np.mean(np.corrcoef(X_new.T))))

    比较结果如下:

    可以看出,经过主成分分析,我们得到了比较好的降维数据,这又一次说明了主成分分析的重要性;

     

    以上就是关于Python和R中主成分分析基础降维功能的介绍,如有不正确之处望指出。

    转载于:https://www.cnblogs.com/feffery/p/8688527.html

    展开全文
  • PCA主成分分析(下)

    2020-08-06 07:07:50
    美,是在高潮处陡然消逝,不落凡尘。数学中的美,是不是也是寻找那个导数为零的极值点?实际问题中,我们认为凸型函数函数中是相对完美而且最容易求极值点的。哦……可惜数学实际上没那么多想象的浪...
  • R-数据挖掘-主成分分析PCA(二) R-数据挖掘-关联规则(三) R-数据挖掘-决策树ID3(四) R-数据挖掘-贝叶斯分类(五) R-数据挖掘-聚类Kmeans(六) R-数据挖掘-聚类DBSCAN(七) 全文逻辑:(读者可将所有...
  • matlab中统计工具箱函数名大全

    千次阅读 2014-11-29 22:18:23
    MATLAB统计工具箱包括概率分布、方差分析、假设检验、分布检验、非参数检验、回归分析、判别分析、主成分分析、因子分析、系统聚类分析、K均值聚类分析、试验设计、决策树、多元方差分析、统计过程控制和统计图形...
  • 首先利用主成分分析技术对功耗数据进行处理,降低数据特征相关性;然后通过遗传算法提高惩罚系数和核函数参数的选择;最后进行硬件木马分类器的构建。实验结果表明,优化SVM方法提高了硬件木马分类器的检测速度和...
  • 目录监督学习三要素模型评估与选择模型介绍感知机k近邻朴素贝叶斯法决策树CART算法逻辑回归支持向量机集成学习adaboost梯度提升树(gbdt)bagging随机森林gbdt于xgboost的区别神经网络多层前馈神经网络非监督学习...
  • R-数据挖掘-主成分分析PCA(二) R-数据挖掘-关联规则(三) R-数据挖掘-决策树ID3(四) R-数据挖掘-贝叶斯分类(五) R-数据挖掘-聚类Kmeans(六) R-数据挖掘-聚类DBSCAN(七) 全文逻辑:(读...
  • 这里有四个比较具体的研究人员试图了解大型数据集的方法以及一些常用的算法。这些方法和算法有许多变数,但是这个表单至少是个良好的...线性投影最流行的做法似乎是主成分分析法。如拓扑学,不过,我们可以想像线...
  • MATLAB统计工具箱

    2015-08-08 10:21:38
    MATLAB统计工具箱包括概率分布、方差分析、假设检验、分布检验、非参数检验、回归分析、判别分析、主成分分析、因子分析、系统聚类分析、K均值聚类分析、试验设计、决策树、多元方差分析、统计过程控制和统计图形...
  • 矩阵的相似标准 线性映射的体积膨胀系数 例子 例题:秩 例子:相似标准 小结 方阵的相合变换 相合不变量 方阵的正交相似变换 方阵的正交相似标准阵 主成分分析 PCA的步骤 长方矩阵的奇异值分解 多元函数的二阶...
  • Penrose伪逆2.10 迹运算2.11 行列式2.12 实例:主成分分析第三章 概率与信息论3.1 为什么要使用概率3.2 随机变量3.3 概率分布3.3.1 离散变量和概率质量函数3.3.2 连续变量和概率密度函数3.4 边缘概
  • 机器学习算法详解

    2018-04-12 15:21:14
    ◦ 4、S型函数(即) ◦ 5、映射为多项式 ◦ 6、使用的优化方法 ◦ 7、运行结果 ◦ 8、使用scikit-learn库中的逻辑回归模型实现 ▪ 逻辑回归_手写数字识别_OneVsAll ◦ 1、随机显示100个数字 ◦ 2、OneVsAll ◦ 3、...
  • R-数据挖掘-主成分分析PCA(二) R-数据挖掘-关联规则(三) R-数据挖掘-决策树ID3(四) R-数据挖掘-贝叶斯分类(五) R-数据挖掘-聚类Kmeans(六) R-数据挖掘-聚类DBSCAN(七) 全文逻辑: ...
  • R-数据挖掘-主成分分析PCA(二) R-数据挖掘-关联规则(三) R-数据挖掘-决策树ID3(四) R-数据挖掘-贝叶斯分类(五) R-数据挖掘-聚类Kmeans(六) R-数据挖掘-聚类DBSCAN(七) 全文逻辑: ...
  • R-数据挖掘-主成分分析PCA(二) R-数据挖掘-关联规则(三) R-数据挖掘-决策树ID3(四) R-数据挖掘-贝叶斯分类(五) R-数据挖掘-聚类Kmeans(六) R-数据挖掘-聚类DBSCAN(七) 全文逻辑: ...
  • R-数据挖掘-主成分分析PCA(二) R-数据挖掘-关联规则(三) R-数据挖掘-决策树ID3(四) R-数据挖掘-贝叶斯分类(五) R-数据挖掘-聚类Kmeans(六) R-数据挖掘-聚类DBSCAN(七) 全文逻辑:(读者可将所有...
  • R-数据挖掘-主成分分析PCA(二) R-数据挖掘-关联规则(三) R-数据挖掘-决策树ID3(四) R-数据挖掘-贝叶斯分类(五) R-数据挖掘-聚类Kmeans(六) R-数据挖掘-聚类DBSCAN(七) 全文逻辑:(读者可将所有...
  • 9.4 主成分分析 9.5 因子分析 9.6 多元方差分析 第10章 统计过程控制 第11章 试验设计 第12章 统计图 第13章 文件输入/输出 第14章 统计演示 第二篇 优化工具箱 第15章 优化工具箱概述 第16章 无约束最优化问题 第17...
  • 现代统计学与SAS应用

    2008-12-01 14:52:34
     第3节 用PRINCOMP过程实现主成分分析  第4节 合成资料的主成分分析 第2章 因子分析  第1节 基本概念  第2节 因子模型  第3节 因子分析的基本定理与任务  第4节 用FACTOR过程实现因子分析 第...
  • 挖掘建模

    2019-04-15 23:43:00
    回归分析是确定预测属性(数值)与其他变量间相互依赖的定量关系最常用的统计学方法。包括线性回归,非线性回归,Logistic回归,岭回归,主成分回归,偏最小二乘回归等模型 决策树 决策树采用自顶向下的递归方式...
  • 14.7 主成分分析法的应用建议 第九周 十五、异常检测(Anomaly Detection) 15.1 问题的动机 15.2 高斯分布 15.3 算法 15.4 开发和评价一个异常检测系统 15.5 异常检测与监督学习对比 15.6 选择特征 15.7 多元高斯...
  • 该方法学方法最初是基于多标准分析(因子分析)来捕获形成城市开放空间劣势和动态的特征,以获得居住满意度的因素。 在对满意度因子(主成分)进行评估的基础上,进行了等级分类,探讨了居民对生物气候更新项目贡献...
  • 3、主成分分析PCA与线性回归的区别 4、PCA降维过程 5、数据恢复 6、主成分个数的选择(即要降的维度) 7、使用建议 8、运行结果 9、使用scikit-learn库中的PCA实现降维 七、异常检测 Anomaly Detection 1、高斯...
  • envi 教程(适用于初学者)

    热门讨论 2009-05-11 11:09:59
    常规处理、几何校正、定标、多光谱分析、高光谱分析、雷达分析、地形地貌分析、矢量应用、神经网络分析、区域分析、GPS联接、正射影象图生成、三维图像生成、丰富的可供二次开发调用的函数库、制图、数据输入/输出等...
  • 2.针对人脸图像在单训练样本下难以被准确识别的问题,提出了一种基于核主成分分析网络(Kerne1 Principle Component Analysis Networks,KPCANet)模型的二阶段投票人脸识别方法。该方法在不使用额外样本数据的情况下,...

空空如也

1 2
收藏数 35
精华内容 14
关键字:

函数型主成分分析