精华内容
下载资源
问答
  • 函数型数据分析的全部代码和数据以及参考文献,全部是本人亲自收集和处理的,共300多M,相关介绍详见我的博文https://blog.csdn.net/lusongno1/article/details/89305520#comments_14878182。
  • 函数型数据主成分分析

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

    函数型数据主成分分析

    基本思想

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

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

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

    故而,求解第一个主成分就变成了求解一个优化问题:
    max ⁡ 1 n ∑ i = 1 n ξ i 2 = max ⁡ 1 n ∑ i = 1 n ( ∫ T β ( s ) x i ( s ) d s ) 2  s.t.  ∥ β ∥ 2 = ∫ T β ( s ) β ( s ) d s = 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} maxn1i=1nξi2 s.t. =maxn1i=1n(Tβ(s)xi(s)ds)2β2=Tβ(s)β(s)ds=1
    求解这个优化问题,我们就得到了第一主成分 β 1 ( s ) \beta^1(s) β1(s)

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

    max ⁡ 1 n ∑ i = 1 n ξ i 2 = max ⁡ 1 n ∑ i = 1 n ( ∫ T β ( s ) x i ( s ) d s ) 2  s.t.  ∥ β ∥ 2 = ∫ T β ( s ) β ( s ) d s = 1 ∫ T β ( s ) β l ( s ) d s = 0 , l = 1 , ⋯   , k − 1 \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} maxn1i=1nξi2=maxn1i=1n(Tβ(s)xi(s)ds)2 s.t. β2=Tβ(s)β(s)ds=1Tβ(s)βl(s)ds=0,l=1,,k1

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

    v ( s , t ) = 1 n − 1 ∑ i = 1 n ( x i ( s ) − x ‾ ( s ) ) ( x i ( 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) v(s,t)=n11i=1n(xi(s)x(s))(xi(t)x(t))

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

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

    定义积分变换:

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

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

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

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

    F V E = ∑ i = 1 K λ i / ∑ i = 1 n − 1 λ i \mathrm{FVE}=\sum_{i=1}^{K} \lambda_{i} / \sum_{i=1}^{n-1} \lambda_{i} FVE=i=1Kλi/i=1n1λi

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

    特征问题求解

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

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

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

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

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

    我们记 X = ( x 1 , x 2 , … , x N ) ′ , Φ = ( Φ 1 , … , Φ k ) ′ , C = ( c i k ) 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=(x1,x2,,xN),Φ=(Φ1,,Φk),C=(cik)N×K

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

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

    定义K阶对称矩阵

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

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

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

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

    将其代入

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

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

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

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

    W W W做cholesky分解,可得 W = L L ′ W=LL' W=LL

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

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

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

    关于基函数的选取

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

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

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

    B j , 0 ( x ) = { 1 , t j ≤ x &lt; t j + 1 0 , else B i , k ( x ) = x − t i t i + k − t i B i , k − 1 ( x ) + t i + k + 1 − x t i + k + 1 − t i + 1 B i + 1 , k − 1 ( x ) , k &gt; 0 \begin{array}{c}{B_{j, 0}(x)=\left\{\begin{array}{l}{1, t_{j} \leq x&lt;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&gt;0}\end{array} Bj,0(x)={1,tjx<tj+10,elseBi,k(x)=ti+ktixtiBi,k1(x)+ti+k+1ti+1ti+k+1xBi+1,k1(x),k>0

    程序代码

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

    主成分分析(PCA)

    数据降维简介

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

    PCA算法的原理解释

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

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

    PCA算法的数学推导

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

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

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

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

    PCA算法简单描述

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

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

    • X − x 0 : = [ x 1 − x 0 , x 2 − x 0 . . . x n − x 0 ] X-x_0 : = [x_1-x_0,x_2-x_0...x_n-x_0] Xx0:=[x1x0,x2x0...xnx0]做奇异值分解: X − x 0 = U Λ V T X-x_0 = U\Lambda V^T Xx0=UΛVT

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

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

    在这里插入图片描述

    PCA算例一

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

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

    X − x 0 = U Λ V T X-x_0 = U\Lambda V^T Xx0=UΛVT

    其中,

    那么,有

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

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

    PCA算例二

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

    其他数据降维方法

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

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

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

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

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

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

    函数型数据主成分分析

    Idea的萌生

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

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

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

    FPCA简介和理论推导

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

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

    • 对函数进行SVD离散化

    • 对函数进行基函数展开

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

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

    常用的基函数有傅里叶基函数和B样条基函数,傅里叶基函数适用于周期性函数数据,B样条基函数适用于非周期函数数据,当然,也可以用多项式基函数。
    B样条基函数的递归定义为:
    B j , 0 ( x ) = { 1 , t j ≤ x < t j + 1 0 , else B i , k ( x ) = x − t i t i + k − t i B i , k − 1 ( x ) + t i + k + 1 − x t i + k + 1 − t i + 1 B i + 1 , k − 1 ( 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} Bj,0(x)={1,tjx<tj+10,elseBi,k(x)=ti+ktixtiBi,k1(x)+ti+k+1ti+1ti+k+1xBi+1,k1(x),k>0
    附录中有一段简单的以多项式为基的MATLAB代码。

    FPCA和PCA的区别和联系

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

    基于FPCA的模型约化

    Onsager原理简介

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

    简单例子:斜板液滴滑动

    问题描述

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

    原理的应用

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

    算法步骤

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

    • 能量耗散函数:
      Φ [ a ˙ , a ] = 1 2 ∫ a 1 a 2 d x ∫ − Y Y d y 3 η h ( V 2 + y 2 W 2 ) \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) Φ[a˙,a]=21a1a2dxYYdyh3η(V2+y2W2)
      其中,
      V ( x , t ) = − 1 H Y ∫ a 1 x ( H ˙ Y + H Y ˙ ) d x W = 1 Y ( Y ˙ + V ∂ x Y ) \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} V(x,t)W=HY1a1x(H˙Y+HY˙)dx=Y1(Y˙+VxY)
      h ( x , y , t ) = H ( x , t ) [ 1 − ( y Y ( x , t ) ) 2 ] h(x, y, t)=H(x, t)\left[1-\left(\frac{y}{Y(x, t)}\right)^{2}\right] h(x,y,t)=H(x,t)[1(Y(x,t)y)2]
      H ( x , t ) = ( x − a 1 ( t ) ) ( a 2 ( t ) − x ) ( a 3 ( t ) + a 4 ( t ) x ) Y ( x , t ) = ( x − a 1 ( t ) ) 1 2 ( a 2 ( t ) − x ) 1 2 ( a 5 ( t ) + a 6 ( 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} H(x,t)=(xa1(t))(a2(t)x)(a3(t)+a4(t)x)Y(x,t)=(xa1(t))21(a2(t)x)21(a5(t)+a6(t)x)

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

    • 势能函数:
      A ( a ) = ∫ a 1 a 2 d x ∫ − Y Y d y [ 1 2 γ θ e 2 + 1 2 γ [ ( ∂ x h ) 2 + ( ∂ y h ) 2 ] + 1 2 ρ g h 2 sin ⁡ α − ρ g x h cos ⁡ α ] \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} A(a)=a1a2dxYYdy[21γθe2+21γ[(xh)2+(yh)2]+21ρgh2sinαρgxhcosα]

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

    数值实验

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

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

    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))
    
    
    展开全文
  • 函数型数据分析若干方法及应用【摘要】:许多领域所收集到的样本观测数据在数据空间中会呈现出一种明显的函数型特征,更多表现形式为光滑的曲线或连续的函数。在对函数型数据进行分析时,将观测到的数据(函数)看作一个...

    函数型数据分析若干方法及应用

    【摘要】:许多领域所收集到的样本观测数据在数据空间中会呈现出一种明显的函数型特征,更多表现形式为光滑的曲线或连续的函数。在对函数型数据进行分析时,将观测到的数据(函数)看作一个整体,而不是一串数字,这是函数型数据分析方法与传统统计分析方法的本质区别。函数型数据分析的目标和传统统计学分析的目标基本一致,即:选用对进一步分析有利的方法来描述数据;为突出不同特征而对数据进行展示;研究数据类型的重要来源和数据之间的变化;利用输入(自变量的信息)来解释输出(因变量)的变化情况;对两组或者更多的某种类型的变量数据进行比较分析等等。

    本文围绕上述目标,构建了一个比较完善、系统的函数型数据分析方法的理论框架。并在函数型数据分析方法的应用方面有所突破和拓展。主要从以下几个方面进行深入地研究。

    1、基函数的选择。基函数展开是将观测得到的离散数据点平滑法预处理的重要手段。我们比较了最常用的傅里叶基函数、Bernstein基函数、B-样条基函数的各自优缺点以及它们的适用范围。更重要的是,本文构造出另外一些基函数:非均匀有理B-样条基函数、三角基函数和混合Bezier类基函数。并且指出,非均匀有理B-样条基函数克服了常用基函数无法准确描述曲线剧烈波动程度的缺点。三角基函数,提供了一个更加广阔的思路。利用该基函数构造的函数型数据曲线在满足一、二阶连续性方面减弱了多项式曲线要求的条件。混合Bezier类基函数不但保留了多项式基函数、有理B-样条基函数的优点,而且使计算更加简单。由此生成的函数型数据曲线具有更大的灵活性,得到更优的结果。

    2、函数型数据曲线的光顺处理。惩罚函数法是传统的函数型数据曲线的光顺方法。除此之外,本文给出了能量光顺法的理论和方法。并在应用B-样条基函数的前提下,深入研究了局部能量最优光顺法、基于曲率均化的B-样条曲线能量光顺法、非均匀样条曲线光顺法以及分层能量光顺法。

    基于均匀节点的B-样条曲线所提出的局部最优光顺法对于局部修改具有重要的实践意义,而且结果表明,该方法计算快速有效。

    基于曲率均化的B-样条曲线能量光顺法,除了考虑数据曲线的应变能满足最小化的要求之外,加入对数据曲线的曲率要求。即对数据曲线施行曲率重新均化,这种均化的方法可以是基于累加弦长的,也可以是基于累加弧长的。

    考虑到在实际数据观测中,所获得的数据未必是均匀的或者等距的。因此,由非均匀节点生成的样条曲线看成是非均匀样条曲线。非均匀样条曲线光顺法就是针对这种情况提出的能量光顺法。

    分层能量法的基本思想是:首先,在一定误差范围内对原始曲线进行小波分解,目的是对曲线进行数据压缩和粗光顺,然后计算边界约束控制点,利用带约束条件的能量法对光顺后的曲线进行边界约束处理和细光顺。

    3、对数据曲线分析方法的研究。在完善了函数型数据的方差分析、典型相关分析之后,主要对函数型数据的主成分分析方法做出了比较全面的研究和实证分析。本文另外一个创新点在于首次推导出基于欧式距离聚类分析的一些结论,指出在使用Fourier基函数和正交基函数展开后,函数型数据曲线的聚类分析可以转化,即是对展开项的生成空间中的点做普通的聚类分析。这大大简化了计算,提高了计算速度。并将Pearson相似系数引入聚类分析之中。特别提出了加权分段聚类分析的概念,不但把握数据曲线簇简单局部形态特征而且提高了聚类质量。

    4、本文从三个角度研究了函数型数据分析方法的应用,这也是对应用领域的一个拓展和突破。首先,分析了汇改前后人民币汇率变动对国内物价的传递效应,将函数型数据分析方法应用到经济领域,并进行了实证分析,得出结论。其次,在基于函数型数据分析方法的笔迹甄别研究中,通过提取数据特征曲线,建立具有良好拟合性质的动态模型,得到了关于笔迹识别一些有用的结论。最后将函数型数据分析方法应用到猪肉价格指数的季节变动中,分析了猪肉价格指数较大波动区间和平稳区间。研究发现,猪肉价格指数的季节变动呈现循环变化的特征。并依据分析结果给出政策建议。

    本文的主要创新包括:

    1、首次将NURBS, Berstein引入到函数想数据分析方法中;

    2、首次构造了三角基函数、混合Bezier类函数并将其应用到函数型数据分析方法中;

    3、首次将能量光顺的思想引入到函数型数据分析方法中,并研究了几种能量光顺的具体实现方法;

    4、指出函数型数据聚类分析方法在正交基下可以通过降维得到实现,并首次提出加权分段聚类的思想;

    5、将函数型数据分析方法应用到经济领域,并进行了实证分析,得出结论。是对应用领域的一个拓展和突破。

    bbb5f0df4a986677d68baa52bb1c42e6.png

    【相似文献】

    中国期刊全文数据库

    前20条

    1

    周妍;;多尺度有限元法的基本原理[J];科协论坛(下半月);2009年04期

    2

    樊晓红;;样条函数与基于样条函数的数据拟合方法[J];科技信息;2011年09期

    4

    马庆华;;某些代数插值的统一与基函数的选取[J];惠州学院学报;1984年S2期

    6

    檀结庆;二次基函数的构造及其应用[J];合肥工业大学学报(自然科学版);1994年02期

    8

    罗军,何国庚,张师帅,黄素逸;高阶平滑基函数重建三维温度场[J];华中理工大学学报;1999年01期

    9

    谭茂金;张庚骥;赵文杰;;普通电阻率测井的新型数值模式匹配解法[J];石油天然气学报;2005年S1期

    10

    李善坡;隋允康;宇慧平;;采用基函数近似的响应面方法[J];北京工业大学学报;2006年S1期

    11

    李小秋;卢俊;贾宏燕;高劲松;冯晓国;孙连春;;具有双频段的十字形复合单元频率选择表面[J];红外与毫米波学报;2007年02期

    12

    陆利正;汪国昭;;二次带形状参数双曲B样条曲线[J];高校应用数学学报A辑;2008年01期

    13

    王阿霞;马逸尘;;复杂起伏地表的POD数值模拟[J];地球科学与环境学报;2008年02期

    14

    王沛;杨嘉林;郭翔;;小波去噪理论的仿真研究[J];科技广场;2009年09期

    15

    余小磊;唐烁;;三角网格上新的插值公式构造[J];中国科学技术大学学报;2011年06期

    16

    宁世光;宣仲良;;量子化学从头计算程序中基函数对称变换表的编制[J];山东师范大学学报(自然科学版);1984年02期

    17

    沈燮昌;Hermite 插值基函数正交展开求和[J];数学学报;1989年01期

    18

    孙红兵,奚梅成;一般的Hermite插值基函数的显式表示[J];中国科学技术大学学报;2001年04期

    19

    龙述尧,陈莘莘;关于无单元法中的插值基函数选取的探讨[J];湖南大学学报(自然科学版);2002年05期

    20

    吉家锋,张晓丽,谢维成;用MATLAB计算等距三次样条插值问题[J];四川工业学院学报;2003年S1期

    中国重要会议论文全文数据库

    前10条

    1

    周志刚;;故障检测中的人工神经网络[A];信号与信息处理技术第三届信号与信息处理全国联合学术会议论文集[C];2004年

    2

    李善坡;隋允康;;基函数近似响应面方法的研究[A];北京力学学会第12届学术年会论文摘要集[C];2006年

    3

    张建伟;玄书鹏;;浸没循环撞击流反应器压力波动信号的Hilbert-Huang谱分析[A];第三届全国化学工程与生物化工年会论文摘要集(上)[C];2006年

    4

    张传东;宋文武;;线天线输入导纳的矩量法分析[A];第六届全国电磁兼容性学术会议2004EMC论文集[C];2004年

    5

    高荣奉;;基础分析的改进有限条法[A];第五届全国结构工程学术会议论文集(第一卷)[C];1996年

    6

    罗健旭;邵惠鹤;;超闭球CMAC学习算法的改进和收敛性分析[A];2004中国控制与决策学术年会论文集[C];2004年

    7

    丁建军;朱剑;陈如山;丁大志;樊振宏;;一种用于电磁散射的多分辩预曲面RWG基函数[A];2009年全国微波毫米波会议论文集(下册)[C];2009年

    8

    杨颖怡;赵延文;陆田;;基于样条时间基函数的精确稳定时间步进算法[A];2009年全国天线年会论文集(上)[C];2009年

    9

    陶诗飞;郭星辰;胡小情;陈如山;;高阶相位基函数结合射线快速多极子方法分析电大尺寸结构问题[A];2011年全国微波毫米波会议论文集(下册)[C];2011年

    10

    杨伍琳;郭建道;黄新媚;罗祖红;;自动气象站数据分析与状态检测报警系统[A];中国气象学会2005年年会论文集[C];2005年

    中国博士学位论文全文数据库

    前10条

    1

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

    2

    殷瑞飞;数据挖掘中的聚类方法及其应用[D];厦门大学;2008年

    4

    徐小红;图像信息的基函数表示方法研究[D];合肥工业大学;2009年

    6

    韩晓刚;复杂运输机械装备运行安全的研究[D];重庆大学;2006年

    7

    崔旭明;B样条与无条件稳定的子区间法[D];天津大学;2007年

    9

    中国硕士学位论文全文数据库

    前10条

    1

    张帆;陀螺加速寿命试验研究[D];南京理工大学;2006年

    3

    朱慧平;智力因素对俄语学习策略影响的实证研究[D];东北师范大学;2009年

    5

    胡勇;基于神经网络的个人信用模型[D];西南财经大学;2006年

    6

    刘利华;信号完整性分析及仿真波形自动分析技术研究[D];西安电子科技大学;2007年

    7

    陈红权;土地利用更新调查技术的研究与应用[D];河海大学;2007年

    8

    郭莉霞;科技型企业员工满意度实证研究[D];北京邮电大学;2008年

    9

    田禾;全国林产品对外贸易信息管理系统的研建[D];北京林业大学;2008年

    10

    范淑静;火焰辐射光谱数据存储、查询及分析系统[D];西安电子科技大学;2008年

    中国重要报纸全文数据库

    前10条

    1

    杨小刚;多一些数据说话,少一些观念争论[N];第一财经日报;2008年

    2

    邹小丹;推出全新商务差旅数据分析[N];中华工商时报;2006年

    3

    本报记者  徐效鸿;指数连阳 板块与个股普涨[N];中国证券报;2006年

    4

    本报记者  徐效鸿;题材股逞强 复牌G股走弱[N];中国证券报;2006年

    5

    ;增长持续快速应用更趋多样[N];人民邮电;2005年

    6

    记者 徐琦;国家环保总局通报松花江水污染情况[N];中国环境报;2005年

    7

    严吕勇;从经验管理到数据管理[N];计算机世界;2005年

    8

    本报记者 徐效鸿;绩优绩差齐跌 板块表现平淡[N];中国证券报;2006年

    9

    惠正一;脑袋VS数据 精准计量营销将效果最大化[N];第一财经日报;2007年

    10

    北京市工商局副局长 罗文阁;深化年检数据分析 拓展企业年检功效[N];中国工商报;2007年

    展开全文
  • 多变量分析中的最大问题莫过于多元线性问题,SPSS降维分析中的主成分分析可以很好地解决这个问题。所谓主成分分析(PCA)也称主分量分析,是有Karl Pearson在1901年提出的,它旨在利用把多个变量指标转化为为少数几...

    810505478457b775f6585ea48e680b08.png

    多变量分析中的最大问题莫过于多元线性问题,SPSS降维分析中的主成分分析可以很好地解决这个问题。所谓主成分分析(PCA)也称主分量分析,是有Karl Pearson在1901年提出的,它旨在利用把多个变量指标转化为为少数几个综合指标,是问题的分析变得更加容易。

    未经许可请勿转载

    更多数据分析内容参看这里

    一. 相关理论

    1. 基本原理

    将多个变量指标通过线性变换浓缩为少数几个主成分指标的多元统计方法。基本思想是把原来多个相关性较强的变量,重新整合为一组互不相关的新的综合指标来代替原来的变量。借助于一个正交变换,将其分量相关的原始随机向量转换成分量不相关的新随机向量。在代数上表现为将原随机向量的协方差阵变换成对角型阵,在几何上表现为将原坐标系变成新的正交坐标系,使之指向样本点散步最开的P个正交方向,然后对多维变量系统进行降维处理,使之能以一个较高的精度转换成低维变量系统,再通过构造适当的价值函数,进一步把低维系统转化为一维系统。方差较大的几个新变量能综合反映原来多个变量包含的主要信息。这几个新变量就是主成分。

    2. 主成分数量筛选依据

    (1)累积方差贡献率:当前m个主成分的累积方差贡献率达到某一特定值(一般80%以上),就可以保留前m个主成分

    (2)特征值:一般选取特征值大于等于1的主成分

    (3)碎石图:一般选取碎石图的曲线上由陡峭变为舒缓的结点前的碎石为主成分

    3. 主成分分析中的主要统计量

    (1)方差贡献率:指的是一个主成分所能够解释的方差占全部方差的比例,这个值越大,说明主成分综合原始变量的信息的能力越强。

    方差贡献率的计算公式为:

    相应的,主成分筛选中所确定的前m个主成分所能解释的全部方差占总方差的比例称为累计方差贡献率。其公式为:

    第一主成分的方差贡献率最大,他能解释原始变量X1,X2....,Xp的能力最强,第2,第3,...第p个主成分的解释能力一次递减。

    (2)特征值:衡量主成分影响力的重要指标,它代表引入该主成分可以解释平均多少原始标量的信息。求出特征值后要按大小予以排列:

    。如果特征值小于1,表示该主成分的解释力非常低,一般以特征值大于1位筛选主成分的标准。

    4. 主成分分析中的适合度检验

    (1)Bartlett球形检验

    原始变量间存在相关性是进行主成分分析的首要条件,否则原始变量无法进行降维处理。为了检验变量之间是否存在相关性,Bartlett在1950年提出了著名的Bartlett球形检验方法,用于检验变量相关系数矩阵是否为单位矩阵。

    设变量相关系数矩阵为R,Bartlett球形检验的统计量为:

    ,其中p为原始变量个数,n 为样本容量,ln|R|是相关系数矩阵行列式的自然对数。

    如果原始变量之间相互独立,那么他们的相关系数矩阵R接近一个单位矩阵。此时|R|接近1,ln|R|约等于0;如果原始变量之间有相关关系,|R|接近于0,ln|R|趋向于负无穷。此时SPSS输出相关系数矩阵的行列式的值,即|R|值。

    Bartlett球形检验的假设是:

    H0:相关系数矩阵是单位矩阵(变量不相关)

    H1:相关系数矩阵不是单位矩阵(变量相关)

    SPSS将输出Bartlett球形检验的卡方统计量,自由度和显著性值。如果显著性值P≤0.05,则认为相关系数矩阵不是单位矩阵,可以进行主成分分析。同时卡方值越大,说明变量之间的相关性越强。

    (2)KMO取样适合度检验统计量

    KMO通过比较样本间的相关系数平方和和偏相关系数平方和的大小以检验赝本是否适合进行主成分分析。如果变量之间的相关系数绝对值较大,而偏相关系数绝对值较小,则表明变量之间高度相关可能与第三变量有关,存在多元线性相关的可能性较大,适合进行主成分分析或因子分析。KMO统计量的定义如下

    , 其中
    为第i个变量和第j个变量的简单相关系数,Pij为第i个变量和第j各个变量在控制了剩余变量后的偏相关系数。

    KMO统计量MSA的取值为0-1,越接近1说明变量间的相关性越强而偏相关性度越低,样本数据越适合做主成分分析和因子分析。根据Kaiser的研究经验,MSA>0.9表示非常合适,0.8-0.9表示合适,0.7-0.8表示一般,0.6-0.7表示尚可,0.5-0.6表示不太合适,0.5以下表示极不合适。

    二. 主成分分析案例应用

    欢迎参加这次双十一狂欢节购物状况调查满意度调查,我们是一个专注于双十一狂欢节购物服务研究的团队。感谢您同意分享您关于双十一狂欢节购物服务的想法。这次调查需要您5-8分钟宝贵的时间。本问卷采用不记名的形式,请您按照实际感受如实回答以下问题。如果您对这次研究的总体报告感兴趣,请您在调查留下电子邮箱地址。谢谢! 本调查包括74个题目,从多方面,多角度对调查主题进行了评测。

    现在以该调查为依据进行主成分分析,尝试从众多变量指标中抽取出能够包含原始便利阿宁大多数信息的少数综合指标,以简化问题的研究。

    1. 导入数据

    c6dd3d5efe0ce3278ca6e66ab0c8e2c2.png

    2. 选择变量

    选择“分析”|“降维”|“因子分析”,打开因子分析对话框。选择左侧框中变量,单击=》箭头,将其选入“变量”列表框。注意主成分分析和因子分析适用的是数值型变量,因此选入分类数据。

    1538f42b92fc9e6ee2aff104f3d848ee.png

    3. 统计量设置

    单击“描述”按钮,打开“因子分析:描述统计”对话框。该对话框包括“统计”和"相关性矩阵”两个选项组。本例子选择默认设置“原始分析结果”和“KMO和Bartlett 球形检验”。单击“继续”按钮返回。

    99ef8a89c4a7b5505952ea3e6d0a507a.png

    4. 主成分抽取设置

    单击“提取”按钮,打开“因子分析:提取”对话框,如下设置,然后返回。

    a5367112fe656f02f5215a9393069e7f.png

    5. 旋转方式设置

    单击“旋转”按钮,打开“因子分析:旋转”对话框,该对话框用于因子分析时的旋转方式设置。因此本例子不做选择,保持默认选择无。

    55f267a74180884f4e11efc7cb303518.png

    6. 得分保存设置

    单击“得分”按钮,打开“因子分析:因子得分”对话框。选中“保存为变量”复选框,将主成分分析得分保存于工作数据中,方法选择“回归”;然后选中“显示因子得分系数矩阵”复选框,单击“继续”返回。

    57e3b3e567f67a52c13606f07ea760e7.png

    7. 选项设置

    单击“选项”按钮,打开“因子分析:选项”对话框,该对话框用于对缺失值的处理和系数显示格式的设置。单击“继续”按钮,返回“因子分析”。

    af3183b1409bfcf58c534ce27d213f29.png

    8。 输出分析结果

    对话框设置完毕,单击“确定”按钮,SPSS输出分析结果

    9. 分析结果解读

    (1)下表是KMO和Bartlett检验结果。从中可以看出,KMO取样适合度检验统计量为0.941的高水平,说明案例中变量之间的信息重叠程度很高,相关度很大。Bartlett检验卡方值16104.441,p值0.000,达到及其显著的水平。两种检验都表明,该例子适合主成分分析。

    6b08874f6304a455287f9e7181192668.png

    (2)公因子方差表示各变量所包含的信息能够被提取的主成分所表示的程度,也称为“共同度”。“初始值”表示每个变量演示信息均为1,即100%。而“提取”栏表示该变量的方差能被主成分所表示的程度,可以看出,变量的方差能被主成分解释在90%以上。

    c89076ed11a2d9187b50db3ed3509ff9.png

    (3)下表给出了提取的主成分的方差总解释量。可以看出,每个主成分能解释的方差的比例不同。第1主成分特征值为61.78,解释原始变量的方差比例为90.853.

    98b12fab2d26e8e34d75faf0d070a38c.png

    (4) 下表给出了2个主成分的成分矩阵,也称为因子载荷,实质是指个主成分和个原始变量的相关系数。通过因子载荷矩阵可以得到原始指标变量的线性组合,如P1=a11*F1+a12*F2其中P1为指标变量1,a11、a12分别为与变量P1在同一行的因子载荷,F1、F2分别为提取的公因子。主成分表达式中相关系数值越大,表示该主成分对原始变量的代表性越大。可以看出,第一主成分y1与个原始变量的关系是2个主成分中较大的,说明它对原始变量的解释量最多。

    b9d59089ddaea057628f567d653bbad6.png

    (5)下表是2个主成分的得分系数。根据主成分得分系数,就可以计算出每个受访者在2个主成分上的得分。假设用X1,X2, ....X74表示所有标准化的原始变量,用F1,F2表示主成分得分。那么就可以得出如下的得分函数 F1=0.019X1+0.019X2+......0.019X74. 简单来说,因子载荷和得分的区别,前者使用主成分来表示原始变量,而后者是用原始变量表示主成分。

    144b868ef8b35868b85443d3401af699.png

    (6)下表是成分得分的协方差矩阵,由于2个主成分互不相关,经过标准化,因此成分得分的协方差矩阵为单位矩阵,即主成分之间的协方差均为0,每个主成分的方差均为1.

    de50d66280d0c6fbe627ac530454b10a.png

    (7)最后我们看一下碎石图。它显示了各个主成分的重要程度,横轴为主成分序号,纵轴为特征值大小。系统选取特征值大于1的前两个成分。

    a69e139ee9f5940c795a49c73709d67f.png

    总之,主成分分析本质上是一种矩阵变换,通过矩阵变换,获得主成分的载荷矩阵,并据此推导出主成分的共同度,方差贡献率等统计量,她并不要求各个主成分具有实际意义。

    b068e6e26a053bf63b3a828ee5be922a.png
    展开全文
  • 可用于 主成分分析、R因子分析、简单相应分析 的R语言函数总结
  • 数据分析-主成分分析流程(R语言)

    千次阅读 2019-07-27 19:19:05
    主成分分析(principal component analysis,PCA)是一种降维技术,把多个变量化为能够反映原始变量大部分信息的少数几个主成分 流程环节为: 1、数据预处理。数值,去缺失值, 2、主成分计算。 3、判断要选择的...
  • 主成分分析详解

    2021-02-09 12:10:11
    多变量大样本无疑会为研究和应用提供了丰富的信息,但也在一定程度上增加了数据采集的工作量,更重要的是在多数情况下,许多变量之间可能存在相关性,从而增加了问题分析的复杂性,同时对分析带来不便。如果分...
  • PCA主成分分析数据进行降维

    千次阅读 2020-02-15 23:42:53
    主成分分析(Principal Component Analysis,PCA), 是一种统计方法。通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量,转换后的这组变量叫主成分。 在实际课题中,为了全面分析问题,往往提出...
  • 目录1 总体主成分1.1 主成分的定义与导出1.2 主成分的性质1.3 从相关矩阵出发求主成分2 样本主成分2.1 从S出发求主成分2.2 从R出发求主成分3 相关的R函数以及实例3.1 `princomp`函数3.2 `summary`函数3.3 `loadings`...
  • 一文看懂主成分分析(PCA)

    千次阅读 2020-12-30 14:26:29
    1 背景2 拆解主成分分析步骤2.1 测试数据2.2 为什么要做主成分分析2.3 step1:数据标准化(中心化)2.4 step2:求相关系数矩阵2.5 step3:计算特征值和特征向量2.6 step4:崖低碎石图和累积贡献图2.7 step5:主成分载荷2.8 ...
  • 今天看了用主成分分析简化数据,就顺便用MNIST数据集做了下实验,想直观地看一下效果,并通过完成这个小demo深入理解下原理。我发现“是什么、能做什么、怎么用、效果是什么、原理是什么、优缺点是什么”这样的思路...
  • 文章目录数据集介绍PCA主成分分析1.基本原理2.代码实现逻辑回归-管道-Pipeline代码模型泛化能力分析 数据集介绍 鸢尾花数据集有三个类别,每个类别有50个样本。其中一个类别与另外两个线性可分,另外两个不能线性可...
  • 微信公众号:幼儿园的学霸 最近在点云处理中,需要获取目标点云的最小包围盒(OBB),在网上看到很多利用PCA求解包围盒的代码,但是代码大多都比较简洁,属于PCA的一个应用,而没有原理的具体描述....主成分分析(Princip.
  • 主成分分析的可视化展示

    千次阅读 2020-08-26 23:08:06
        主成分分析法(PCA)就是一种常见的降维算法,它能够降低数据维度,减少高维数据分析难度。主成分分析法能够在实现降维的同时,能够尽量的保证信息损失。因此在很多分析工作中,可以通过提炼主成分的方式,仅...
  • 二维数据主成分分析的原理—— PCA投影PCA与线性回归的区别matlab自带PCA函数代码 数据降维 二维数据->一维数据、 降维到一维之后的数据: 三维数据->二维数据 上图中的数据可以投影到二维平面上,大致...
  • 目录数据集介绍PCA主成分分析1.基本原理2.代码实现逻辑回归-管道-Pipeline代码模型泛化能力分析 数据集介绍 鸢尾花数据集有三个类别,每个类别有50个样本。其中一个类别与另外两个线性可分,另外两个不能线性可分。 ...
  • Matlab程序实现心脏形态学形状主成分的聚类分析.pdf