精华内容
下载资源
问答
  • 利用半方差分析方法,进行区域变量的空间自相关性分析,了解其空间变异特征
  • 均值-CVaR和均值-下半方差投资组合优化研究,舒燕菲,张鹏,对于投资组合风险风险度量,不同学者有不同主张。考虑在市场摩擦条件下,本文选择了下半方差和CVaR度量投资组合的风险,分别提出了�
  • 首先在标准差和下行标准差的基础上计算了贝塔和下行贝塔值,...由于数据相对集中,又利用了Spearman秩相关系数进行了验证,结果表明:在我国这样一个缺少对冲工具的新兴金融市场-k,应用半方差来刻画风险要优于方差。
  • 分为三个函数,分别进行半方差的拟合和克里金插值
  • 半方差分析ppt

    2013-04-14 10:08:27
    系统的介绍了半方差的理论,使初学者能够初步了解半方差分析。
  • 论文研究-基于滚动经济回撤约束和下半方差的最优投资组合策略.pdf, 最大回撤和下半方差均是投资者在实际投资过程中非常关注的风险指标.本文创造性地将这两个风险指标...
  • matlab开发-分形时间周期分析的半方差方法和比例窗方差方法。分形时间序列分析的半方差法和比例窗方差法。
  • 利用半方差分析方法,进行区域变量的空间自相关性分析,了解其空间变异特征
  • 基于MATLAB的变异函数计算与经验半方差图绘制1 数据处理1.1 数据读取1.2 异常数据剔除1.3 正态分布检验及转换2 距离量算3 距离分组4 平均距离、半方差计算及其绘图5 绘图结果   在前期的博客...


      在前期的 博客(https://blog.csdn.net/zhebushibiaoshifu/article/details/113943720)中,我们详细介绍了地学计算的几个基本概念,并对其数学推导公式加以了梳理。接下来,我将通过几篇新的专题博客,对地学计算相关的代码、操作加以实践与详细讲解。本篇博客便是第一篇—— 基于MATLAB的空间数据变异函数计算与经验半方差图绘制
      另一方面,由于上述博客所涉及的相关理论概念较为抽象,往往需要结合实践才可以更好理解,因此大家可以将上述博客与本篇及后期的其它地学计算博客一同来看,可以更好理解相关理论的含义。
      其中,由于本文所用的数据并不是我的,因此遗憾不能将数据一并展示给大家;但是依据本篇博客的思想与对代码的详细解释,大家用自己的数据,可以将空间数据变异函数计算与经验半方差图绘制的全部过程与分析方法加以完整重现。

    1 数据处理

    1.1 数据读取

      本文中,我的初始数据为某区域658个土壤采样点的空间位置(X与Y,单位为米)、pH值有机质含量全氮含量。这些数据均存储于“data.xls”文件中;而后期操作多于MATLAB软件中进行。因此,首先需将源数据选择性地导入MATLAB软件中。
      利用MATLAB软件中xlsread函数可以实现这一功能。具体代码附于“1.3 正态分布检验及转换”处。

    1.2 异常数据剔除

      得到的采样点数据由于采样记录、实验室测试等过程,可能具有一定误差,从而出现个别异常值。选用“平均值加标准差法”对这些异常数据加以筛选、剔除。
      分别利用“平均值加标准差法”中“2S”与“3S”方法加以处理,发现“2S”方法处理效果相对后者较好,故后续实验取“2S”方法处理结果继续进行。
      其中,“2S”方法是指将数值大于或小于平均值±2倍标准差的部分视作异常值,“3S”方法则是指将数值大于或小于平均值±3倍标准差的部分视作异常值。
      得到异常值后,将其从658个个采样点中剔除;剩余的采样点数据继续后续操作。
      本部分具体代码附于“1.3 正态分布检验及转换”处。

    1.3 正态分布检验及转换

      计算变异函数需建立在初始数据符合正态分布的假设之上;而采样点数据并不一定符合正态分布。因此,我们需要对原始数据加以正态分布检验。
      一般地,正态分布检验可以通过数值检验直方图QQ图等图像加以直观判断。本文综合采取以上两种数值、图像检验方法,共同判断正态分布特性。
      针对数值检验方法,我在一开始准备选择采用Kolmogorov-Smirnov检验方法;但由于了解到,这一方法仅仅适用于标准正态检验,因此随后改用Lilliefors检验
      Kolmogorov-Smirnov检验通过样本的经验分布函数与给定分布函数的比较,推断该样本是否来自给定分布函数的总体;当其用于正态性检验时只能做标准正态检验。
      Lilliefors检验则将上述Kolmogorov-Smirnov检验改进,其可用于一般的正态分布检验。
      QQ图(Quantile Quantile Plot)是一种散点图,其横坐标表示某一样本数据的分位数,纵坐标则表示另一样本数据的分位数;横坐标与纵坐标组成的散点图代表同一个累计概率所对应的分位数。因此,QQ图具有这样的特点:
      y=x
      针对这一直线,若散点图中各点均在直线附近分布,则说明两个样本为同等分布;因此,若将横坐标(纵坐标)表示为一个标准正态分布样本的分位数,则散点图中各点均在上述直线附近分布可以说明,纵坐标(横坐标)表示的样本符合或基本近似符合正态分布。本文采用将横坐标表示为正态分布的方式。
      此外,PP图(Probability Probability Plot)同样可以用于正态分布的检验。PP图横坐标表示某一样本数据的累积概率,纵坐标则表示另一样本数据的累积概率;其根据变量的累积概率对应于所指定的理论分布累积概率并绘制的散点图,用于直观地检测样本数据是否符合某一概率分布。和QQ图类似,如果被检验的数据符合所指定的分布,则其各点均在上述直线附近分布。若将横坐标(纵坐标)表示为一个标准正态分布样本的分位数,则散点图中各点均在直线附近分布可以说明,纵坐标(横坐标)表示的样本符合或基本近似符合正态分布。
      三种土壤属性,我选择首先以pH数值为例进行操作。通过上述数值检验、图像检验方法,检验得到剔除异常值后的原始pH数值数据并不符合正态分布这一结论。因此,尝试对原数据加以对数开平方等转换处理;随后发现,原始pH值开平方数据的正态分布特征虽然依旧无法通过较为严格的Lilliefors检验,但其直方图、QQ图的图像检验结果较为接近正态分布,并较之前二者更加明显。故后续取开平方处理结果继续进行。
      值得一提的是,本文后半部分得到pH值开平方数据的实验变异函数及其散点图后,在对其余两种空间属性数据(即有机质含量与全氮含量)进行同样的操作时,发现全氮含量数据在经过“2S”方法剔除异常值后,其原始形式的数据是可以通过Lilliefors检验的,且其直方图、QQ图分布特点十分接近正态分布。
      我亦准备尝试对空间属性数据进行反正弦转换。但随后发现,已有三种属性数值的原始数据并不严格分布在-1至1的区间内,因此并未对其进行反正弦方式的转换。
      经过上述检验、转换处理过后的图像检验结果如下所示。
    在这里插入图片描述
      以上部分代码如下:

    clc;clear;
    info=xlsread('data.xls');
    oPH=info(:,3);
    oOM=info(:,4);
    oTN=info(:,5);
     
    mPH=mean(oPH);
    sPH=std(oPH);
    num2=find(oPH>(mPH+2*sPH)|oPH<(mPH-2*sPH));
    num3=find(oPH>(mPH+3*sPH)|oPH<(mPH-3*sPH));
    PH=oPH;
    for i=1:length(num2)
        n=num2(i,1);
        PH(n,:)=[0];
    end
    PH(all(PH==0,2),:)=[];
     
    %KSTest(PH,0.05)
    H1=lillietest(PH);
     
    for i=1:length(PH)
        lPH(i,:)=log(PH(i,:));
    end
     
    H2=lillietest(lPH);
     
    for i=1:length(PH)
        sqPH(i,:)=(PH(i,:))^0.5;
    end
     
    H3=lillietest(sqPH);
     
    % for i=1:length(PH)
    %     arcPH(i,:)=asin(PH(i,:));
    % end
    % 
    % H4=lillietest(arcPH);
     
    subplot(2,3,1),histogram(PH),title("Distribution Histogram of pH");
    subplot(2,3,2),histogram(lPH),title("Distribution Histogram of Natural Logarithm of pH");
    subplot(2,3,3),histogram(sqPH),title("Distributio n Histogram of Square Root of pH");
    subplot(2,3,4),qqplot(PH),title("Quantile Quantile Plot of pH");
    subplot(2,3,5),qqplot(lPH),title("Quantile Quantile Plot of Natural Logarithm of pH");
    subplot(2,3,6),qqplot(sqPH),title("Quantile Quantile Plot of Square Root of pH");
    

    2 距离量算

      接下来,需要对筛选出的采样点相互之间的距离加以量算。这是一个复杂的过程,需要借助循环语句。
      本部分具体代码如下。

    poX=info(:,1);
    poY=info(:,2);
    dis=zeros(length(PH),length(PH));
    for i=1:length(PH)
        for j=i+1:length(PH)
            dis(i,j)=sqrt((poX(i,1)-poX(j,1))^2+(poY(i,1)-poY(j,1))^2);
        end
    end
    

    3 距离分组

      计算得到全部采样点相互之间的距离后,我们需要依据一定的范围划定原则,对距离数值加以分组。
      距离分组首先需要确定步长。经过实验发现,若将步长选取过大会导致得到的散点图精度较低,而若步长选取过小则可能会使得每组点对总数量较少。因此,这里取步长为500米;其次确定最大滞后距,这里以全部采样点间最大距离的一半为其值。随后计算各组对应的滞后级别、各组上下界范围等。
      本部分具体代码附于本文“4 平均距离、半方差计算及其绘图”处。

    4 平均距离、半方差计算及其绘图

      分别计算各个组内对应的点对个数、点对间距离总和以及点对间属性值差值总和等。随后,依据上述参数,最终求出点对间距离平均值以及点对间属性值差值平均值。
      依据各组对应点对间距离平均值为横轴,各组对应点对间属性值差值平均值为纵轴,绘制出经验半方差图。
      本部分及上述部分具体代码如下。

    madi=max(max(dis));
    midi=min(min(dis(dis>0)));
    radi=madi-midi;
    ste=500;
    clnu=floor((madi/2)/ste)+1;
    ponu=zeros(clnu,1);
    todi=ponu;
    todiav=todi;
    diff=ponu;
    diffav=diff;
    for k=1:clnu
        midite=ste*(k-1);
        madite=ste*k;
        for i=1:length(sqPH)
            for j=i+1:length(sqPH)
                if dis(i,j)>midite && dis(i,j)<=madite
                    ponu(k,1)=ponu(k,1)+1;
                    todi(k,1)=todi(k,1)+dis(i,j); diff(k,1)=diff(k,1)+(sqPH(i)-sqPH(j))^2;
                end
            end
        end
        todiav(k,1)=todi(k,1)/ponu(k,1);
        diffav(k,1)=diff(k,1)/ponu(k,1)/2;
    end
    plot(todiav(:,1),diffav(:,1)),title("Empirical Semivariogram of Square Root of pH");
    xlabel("Separation Distance (Metre)"),ylabel("Standardized Semivariance");
    

    5 绘图结果

      通过上述过程,得到pH值开平方后的实验变异函数折线图及散点图。
    在这里插入图片描述
    在这里插入图片描述
      可以看到,pH值开平方后的实验变异函数较符合于有基台值的球状模型或指数模型。函数数值在距离为0至8000米区间内快速上升,在距离为8000米后数值上升放缓,变程为25000米左右;即其“先快速上升,再增速减缓,后趋于平稳”的图像整体趋势较为明显。但其数值整体表现较低——块金常数为0.004左右,而基台值仅为0.013左右。为验证数值正确性,同样对有机质、全氮进行上述全程操作。
      得到二者对应变异函数折线图与散点图。

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

      由以上三组、共计六幅的pH值开平方、有机质与全氮对应的实验变异函数折线图与散点图可知,不同数值对应实验变异函数数值的数量级亦会有所不同;但其整体“先快速上升,再增速减缓,后趋于平稳”的图像整体趋势是十分一致的。
      此外,如上文所提到的,针对三种空间属性数据(pH值、有机质含量与全氮含量)中最符合正态分布,亦是三种属性数据各三种(原始值、取对数与开平方)、共九种数据状态中唯一一个通过Lilliefors正态分布检验的数值——全氮含量经过异常值剔除后的原始值,将其正态分布的图像检验结果特展示如下。
    在这里插入图片描述
    至此,我们就完成了全部的操作、分析过程~

    欢迎关注公众号:疯狂学习GIS
    在这里插入图片描述

    展开全文
  • 半方差分析,半方差分析GS,matlab源码.zip
  • matlab计算经验半方差(变异函数)

    千次阅读 多人点赞 2020-04-06 11:02:45
    2 其中 hhh为距离(滞后距) N(h)N(h)N(h)为样点中符合该距离 hhh的点对数量 xix_ixi​和 xi+hx_i+hxi​+h为各点对位置 2随机样点的试验半方差 在很多情况下,样点并不是规则采集(机械采样)的,而是呈现不规则分布,...

    前言

    前段时间地统计的课设代码的确有点折磨人,但是还好算法过程比较简单,没有涉及到特别复杂的步骤,matlab在矩阵运算方面有大量的库函数以及其在索引方面的优势,使得我们真正能够实现向量化(矩阵化编程),希望大家看到这篇文章过后能够真正取思考,以前的哪些代码可以将for循环等若干繁琐步骤去掉,减少代码行数,掌握向量化编程的精髓。

    1 变异函数

    这里引用一下杨老师ppt上的一个简单定义,一维条件下,当空间点 x x x在一维 x x x轴上变化时,区域化变量 Z ( x ) Z(x) Z(x)在点 x x x x + h x+h x+h处的值 Z ( x ) Z(x) Z(x) Z ( x + h ) Z(x+h) Z(x+h)差的方差一半定义为区域化变量Z(x)在x轴方向上的变异函数。
    通常区域化变量 Z ( x ) Z(x) Z(x)所描述的现象是二维和三维的,则变异函数可定义为在任一方向 α α α,相距 ∣ h ∣ |h| h的两个区域化变量值 Z ( x ) Z(x) Z(x) Z ( x + h ) Z(x+h) Z(x+h)的增量的方差。
    2 γ ( x , h ) = V a r [ Z ( x ) − Z ( x + h ) ] = E [ ( Z ( x ) − Z ( x + h ) ) 2 ] − ( E [ ( Z ( x ) − Z ( x + h ) ) ] ) 2 2\gamma(x,h)=Var[Z(x)-Z(x+h)]\newline =E[(Z(x)-Z(x+h))^2]-(E[(Z(x)-Z(x+h))])^2 2γ(x,h)=Var[Z(x)Z(x+h)]=E[(Z(x)Z(x+h))2](E[(Z(x)Z(x+h))])2
    其中 γ ( x , h ) γ(x,h) γ(x,h)称为半变异函数,但有时为方便也称变异函数
    变异函数计算公式
    γ ( h ) = 1 N ( h ) ∑ i = 1 n [ Z ( x i ) − Z ( x i + h ) ] 2 \gamma(h)=\frac{1}{N(h)}\sum_{i=1}^n[Z(x_i)-Z(x_i+h)]^2 γ(h)=N(h)1i=1n[Z(xi)Z(xi+h)]2
    其中 h h h为距离(滞后距)
    N ( h ) N(h) N(h)为样点中符合该距离 h h h的点对数量
    x i x_i xi x i + h x_i+h xi+h为各点对位置

    2随机样点的试验半方差

    在很多情况下,样点并不是规则采集(机械采样)的,而是呈现不规则分布,这时,如何计算试验半方差?
    一般地,在实际计算时,假设步长为 l a g lag lag,当前滞后级别为 n n n( n n n为正整数),则 h = n ∗ l a g h=n*lag h=nlag,应该这样处理:

    计算步骤
    1.研究区所有点,找到点对 ( P i , P j ) (P_i,P_j) (Pi,Pj),其符合条件: ( n − 1 ) ∗ l a g < d i s ( P i , P j ) < = n ∗ l a g (n-1)*lag<dis (P_i,P_j)<=n*lag (n1)lag<dis(Pi,Pj)<=nlag,它们之间的距离记为 D I S i DIS_i DISi
    2.计算 [ z ( p i ) − z ( p j ) ] 2 [z(p_i)-z(p_j)]^2 [z(pi)z(pj)]2,记为 S i S_i Si.
    3.设找到 N ( h ) N(h) N(h)个这样的点对,计算平均距离 h a v g = 1 N ( h ) ∑ i = 1 N ( h ) D I S i h_{avg}=\frac{1}{N(h)}\sum_{i=1}^{N(h)}DIS_i havg=N(h)1i=1N(h)DISi
    4.计算 r ∗ ( h a v g ) = 1 2 N ( h ) ∑ i = 1 N ( h ) S i r^*(h_{avg})=\frac{1}{2N(h)}\sum_{i=1}^{N(h)}S_i r(havg)=2N(h)1i=1N(h)Si n n n滞后级别上的经验半方差值。
    5.将各个级别的 ( h a v g , r ∗ ( h a v g ) ) (h_{avg},r*_{(havg)}) (havgr(havg)),绘制在图上,形成经验半方差图

    3代码实例

    3.1数据集

    本次数据来源于杨老师分享的沙洋县土壤样点数据集,如下图
    在这里插入图片描述
    数据集前两列为 ( x , y ) (x,y) (x,y)坐标,后面几列为属性数据

    3.2 代码实例

    这里以全氮数据为例计算经验半方差
    由于数据为EXCEL形式,因此可调用xlsread函数读取文件将数据文件存储为矩阵,而将表中的文本类型忽略
    数据预处理

    data=xlsread('F:\geostatistical_test\test\data_and_materials\data.xls');
    N=data(:,5);
    [h,p,ks]=lillietest(N);%h==0,通过正态分布检验,p>0.05接受检验
    figure;
    histogram(N,10);
    figure;
    qqplot(N);
    aveN=mean(N);
    varN=var(N,1);%除以n的方差
    stdN=sqrt(varN);
    us=aveN+2*stdN;
    ux=aveN-2*stdN;
    logi=(N>us)+(N<ux);
    logi=logical(logi);
    N(logi)=[];
    data(logi,:)=[];
    [h1,p1,ks1]=lillietest(N);%h==0,通过正态分布检验,p>0.05接受检验
    figure;
    qqplot(N);
    
    

    原始数据
    原始数据QQ图
    去除异常值后
    在这里插入图片描述

    %计算出距离矩阵
    D=pdist(data(:,1:2));
    D1=squareform(D);%计算出距离矩阵
    lag=500;
    maxdist=max(max(D1));%计算最大距离
    D2=triu(D1);%将距离矩阵转化成上三角,下三角为0,方便后面符合序号对提取
    %下三角部分全部变为0
    max_lagcount=round(maxdist/(2*lag))+1;%计算出最多可以求得的延迟阶数
    
    

    将处理好的数据计算各个点对所对应的距离矩阵,其中D1就是一个 n ∗ n n*n nn的距离矩阵,在对角线上为0;以500米为步长,找到最大距离(这里用max(max(D1)))是因为max函数是默认对矩阵每一列求最大值,通过计算可以得到最大能计算的步长max_lagcount
    开始进行for循环计算组合的点对

    循环代码如下:
    result=[];%产生一个空矩阵,存储结果点对
    for n=1:max_lagcount
        lower=(n-1)*lag;
        upper=n*lag;
        logi=logical((D2>lower)+(D2<=upper)-1);%找到符合条件的点对逻辑矩阵
        [tag1,tag2]=find(logi);%[tag1,tag2]表示点对序号
        N1=N(tag1);
        N2=N(tag2);
        S=sum((N1-N2).^2);
        Nh=length(tag1);
    	rh=S/(2*Nh);
    	havg=sum(D1((tag2-1).*size(D2,1)+tag1))/Nh;
        %havg=trace(D1(tag1,tag2))/Nh;
        temp=[havg,rh];
        result=[result;temp];
    end
    
    

    其中lower,upper分别为距离限制条件,
    其中logi=logical((D2>lower)+(D2<=upper)-1);%找到符合条件的点对逻辑矩阵,用的是D2上三角矩阵,若用原始距离矩阵(是对称矩阵),会将点对所有点对重复记录一次,因此选取上三角部分就可记录一次,因为下三角全为0

    这里的tag1,tag2均为一列矩阵,其中 ( t a g 1 i , t a g 2 i ) (tag1_i,tag2_i ) (tag1i,tag2i)就表示着对应的点对序号,因此可以直接通过序号来对N矩阵索引,求得 ( N 1 − N 2 ) . 2 (N1-N2).^2 (N1N2).2,其中每个值就对应的是 S i S_i Si,然后后续调用sum函数求和,其中Nh=length(tag1),找到目标点对得数量,只需要调用length求矩阵长度得函数即可,因此求得 r ∗ ( h a v g ) r^* (h_{avg} ) r(havg)对应于rh
    havg=sum(D1((tag2-1).*size(D2,1)+tag1))/Nh;这里为什么用
    D1((tag2-1).*size(D2,1)+tag1)这么复杂的索引方式对D1距离矩阵进行索引呢,因为前面的 ( t a g 1 i , t a g 2 i ) (tag1_i,tag2_i ) (tag1i,tag2i)表示的是D1中第 t a g 1 i tag1_i tag1i行第 t a g 2 i tag2_i tag2i列的数据,但是如果用D1(tag1,tag2))索引方式的话会将其中每个元素来进行匹配,因此我们要将点对表示的索引转化为一个值即序号索引,因为matlab里面默认情况下是先对第一列索引完后再对第二列索引,
    因此第一列第 i i i行对应于序号为 i i i
    第二列第 i i i行对应于 m + i m+i m+i
    因此第 i i i行第 j j j列对应于序号 ( j − 1 ) ∗ m + i (j-1)*m+i (j1)m+i
    讲的比较通俗一点就是 a ( i , j ) = = a ( ( j − 1 ) ∗ m + i ) a(i,j)==a((j-1)*m+i) a(i,j)==a((j1)m+i)
    m m m为矩阵 a a a的行数,这里举一个简单的例子说明
    在这里插入图片描述
    如果硬要用D1(tag1,tag2)索引的方式也是可行的,
    ,其中D1(tag1,tag2),其中tag1,tag2是目标点对序号得矩阵,但是matlab里面两个索引均为矢量时,会将矢量中每个元素都进行匹配来进行索引,进而返回索引得到得值,这里举一个简单例子说明:
    a=[1,2,3;4,5,6;7,8,9]
    a([1;2],[2;3])
    返回结果分别为
    在这里插入图片描述
    但是我们只想要(1,2),(2,3)(第一行第二列,第二行第三列)处的值,发现其对应的值均会在对角线上,因此对角线上就是序号一一对应处的矩阵索引值,因此这里将调用矩阵的迹trace函数,即对角线元素之和,将这些点对对应的距离相加,最后得到 h a v g h_{avg} havg
    现在问题来了,哪种索引方式更快消耗的内存更小,因此我在代码快的前后加上了计时函数tic,toc根据最终的运行效果对比,基于
    havg=sum(D1((tag2-1).*size(D2,1)+tag1))/Nh
    索引方式所消耗计算时间为
    1.1038s
    基于trace点对索引的方式消耗计算时间为
    2.9563s

    最后将这些点对进行记录
    temp=[havg,rh];result=[result;temp];
    最终会返回一个两列的result矩阵,第一列为 h a v g h_{avg} havg,第二列为 r ∗ ( h a v g ) r^* (h_{avg}) r(havg)
    在这里并没有考虑 N ( h ) = 0 N(h)=0 N(h)=0的情况,假设N(h)=0,那么对应的 ∑ i = 1 N ( h ) D I S i = 0 \sum_{i=1}^{N(h)}DIS_i=0 i=1N(h)DISi=0 ∑ i = 1 N ( h ) S i = 0 \sum_{i=1}^{N(h)}S_i=0 i=1N(h)Si=0,因此 r ∗ ( h a v g ) r^* (h_{avg}) r(havg) h a v g h_{avg} havg均为0/0,在matlab里面计算得到的结果为NaN
    结果为空
    因此我们只需将最终得到的矩阵中将其中为NaN的一行去掉即可
    调用isnan函数,
    result(isnan(result))=[];
    绘图表示:
    plot(result(:,1),result(:,2),’*r’)
    title(‘变异函数散点图’);
    xlabel(‘平均距离’);
    ylabel(‘滞后级别半方差值’);
    在这里插入图片描述
    完整代码如下:

    tic;
    data=xlsread('F:\geostatistical_test\test\data_and_materials\data.xls');
    N=data(:,5);
    [h,p,ks]=lillietest(N);%h==0,通过正态分布检验,p>0.05接受检验
    figure;
    histogram(N,10);
    figure;
    qqplot(N);
    aveN=mean(N);
    varN=var(N,1);%除以n的方差
    stdN=sqrt(varN);
    us=aveN+2*stdN;
    ux=aveN-2*stdN;
    logi=(N>us)+(N<ux);
    logi=logical(logi);
    N(logi)=[];
    data(logi,:)=[];
    [h1,p1,ks1]=lillietest(N);%h==0,通过正态分布检验,p>0.05接受检验
    figure;
    qqplot(N);
    %计算出距离矩阵
    D=pdist(data(:,1:2));
    D1=squareform(D);%计算出距离矩阵
    lag=500;
    maxdist=max(max(D1));%计算最大距离
    D2=triu(D1);%将距离矩阵转化成上三角,下三角为0,方便后面符合序号对提取
    %下三角部分全部变为0
    max_lagcount=round(maxdist/(2*lag))+1;%计算出最多可以求得的延迟阶数
    result=[];
    for n=1:max_lagcount
        lower=(n-1)*lag;
        upper=n*lag;
        logi=logical((D2>lower)+(D2<=upper)-1);%找到符合条件的点对逻辑
        [tag1,tag2]=find(logi);%[tag1,tag2]表示点对序号
        N1=N(tag1);
        N2=N(tag2);
        S=sum((N1-N2).^2);
        Nh=length(tag1);
        rh=S/(2*Nh);
        %havg=trace(D1(tag1,tag2))/Nh;
        havg=sum(D1((tag2-1).*size(D2,1)+tag1))/Nh;
        temp=[havg,rh];
        result=[result;temp];
    end
    result(isnan(result))=[];
    figure;
    plot(result(:,1),result(:,2),'*r')
    title('变异函数散点图');
    xlabel('平均距离');
    ylabel('滞后级别半方差值');
    t1=toc
    

    4总结

    matlab的向量化编程是非常方便的,因为matlab针对矩阵索引支持逻辑索引以及位置索引,在筛出异常值与对特定条件的值进行索引时是非常有效的,其次调用matlab的内置函数一般计算速度要比自己写for循环的最终代码要快得多,如果不相信的话,同学们可以把矩阵化编程的那些语句以及函数改写成for循环试一试比较最终的计算时间,希望大家在今后矩阵运算的编程过程中尽量少用for循环,尽量理解如下代码的含义,方便今后套用:

    logi=(N>us)+(N<ux);
    logi=logical(logi);
    N(logi)=[];
    data(logi,:)=[];
    
    logi=logical((D2>lower)+(D2<=upper)-1);%找到符合条件的点对逻辑
    [tag1,tag2]=find(logi);%[tag1,tag2]表示点对序号
    
    havg=trace(D1(tag1,tag2))/Nh;
    havg=sum(D1((tag2-1).*size(D2,1)+tag1))/Nh;
    

    本次代码全为潘老师自主编写,若要引用,注意引用规范

    撰写本次博客也掌握了对markdown编写数学公式的技巧,也是对自己一个方面的提升,在地统计实验课程结束后会将遗传算法以及克里格插值部分的源代码分享出来造福大家!敬请期待!

    展开全文
  • IDL code: ; Define a procedure to return γ(h) and the partial derivatives according Equation(7) pro Semi_Var,X,A,F,pder F=A[0]+A[1](1.5(X/A[2])-0.5*(X/A[2])^3) if n_params() ge 4 then begin ...

    IDL code:

    ; Define a procedure to return γ(h) and the partial derivatives according Equation(7)
    pro Semi_Var,X,A,F,pder
    F=A[0]+A[1](1.5(X/A[2])-0.5*(X/A[2])^3)
    if n_params() ge 4 then begin
    pder=fltarr(n_elements(X),3)
    pder[,0]=1
    pder[
    ,1]=(1.5*(X/A[2])-0.5*(X/(A[2]^3)
    pder[,2]=1.5A[1]*(X3/(A[2]4)-(X/A[2])^2))
    endif
    end

    ; Compute the fit to the function
    ; Input the independent h and dependent variables r calculated according Equation (6)
    ; Define a vector of weights.
    weights = 1.0/r

    ; Provide an initial guess of the function’s parameters.
    A = [1.0, 5.0, 30.0]

    ; Compute the parameters.
    yfit = CURVEFIT(h, r, weights, A, SIGMA, FUNCTION_NAME=‘Semi_Var’)

    C0=A[0] & C=A[1] & a=A[2]

    Fitting results:
    Location: 29.8065N,121.7873E
    Landsat Data: LT05_L1TP_20091021 (Red Band)
    1.0 Km: C0=3.70000, C=8.44564, a=739.253
    1.5 Km: C0=3.26735, C=5.64539, a=870.949
    2.0 Km: C0=3.33918, C=4.27791, a= 985.988

    展开全文
  • 在考虑系统旋转备用的容量成本和能量成本、因购买旋转备用而减少的停电损失,以及旋转备用效益的离散程度的基础上,引入风险偏好系数,并采用证券投资组合中的加权半方差度量风险。以期望旋转备用效益最大和风险最小...
  • 半方差分析的原理和应用 学习《景观生态学——格局、过程、尺度与等级》(第二版)邬建国著有感 页码范围:P129-136 原理 空间自相关可以用协方差与方差之比表示: 方差 协方差即两个样本间的差异 协方差的实际...

    半方差(semi-variogram)分析的原理和应用

    学习《景观生态学——格局、过程、尺度与等级》(第二版)邬建国著有感
    页码范围:P129-136

    半方差(semi-variogram)又称变异函数!

    1. 原理

    空间自相关可以用协方差与方差之比表示:r(h) = C(h)/C(0)
    协方差即两个样本间的差异在这里插入图片描述
    协方差的实际计算公式:在这里插入图片描述

    半方差的定义在这里插入图片描述
    在这里插入图片描述

    2. 应用【结合《空间数据分析教程》王劲峰】

    以间隔距离h为横坐标,以半方差为纵坐标,得到半方差图。

    • 如何根据半方差图来解释空间结构特征?
      通过分析变化曲线的形状,以此选定用什么样的模型来描述变化。

    我的理解:
    半方差可以用于空间连续数据选定研究单元大小【在景观生态学中就是“斑块大小”】,基于不同单元大小反映的半方差变化,选定后续需要的模型。

    半变异函数在度量空间自相关性时,反映的是全局的空间自相关程度。

    • ArcGIS实现

    Geostatistical Analyst工具条→【Explore Data】→【Semivariogram/Covariance Cloud】
    得到半变异函数图,解读函数图,判断空间自相关显著否。若不显著就不应当使用kriging插值及基于空间自相关性的各种统计。

    展开全文
  • 该算法使用半方差技术计算分形维数。 它有助于评估图像中纹理图案的方向性。 沿水平方向的半方差被定义为在所有像素 N 上的总和,即由 x 方向上的位移 'h' 分隔的像素的强度 I(x,y) 之和。 与垂直方向类似,它是在 y...
  • Matlab半方差函数/变异函数

    千次阅读 2020-05-04 17:33:46
    参考代码:https://blog.csdn.net/qq_44589327/article/details/105338323 原文博主对变异函数的原理和代码部分解释的相当详细,包括正态分布统计、QQplot图、剔除残差后的QQplot图和变异函数曲线。 ...
  • 用于分形时间序列分析的半方差方法和缩放窗口方差方法。 参考:Evaluating scaled windowed variance methods for estimating the Hurst coefficient of time series Physica A: Statistical and Theoretical ...
  • 1.半方差函数:也称空间变异函数是地统计学的重要组成部分,是抽样间隔为h时样本值方差数学期望的一半。以变异函数K(h)为Y轴,抽样间隔h为x轴,可绘成变异函数曲线图。 2.块金值(Nugget)用Co表示:也叫块金方差...
  • GS+7.0半方差函数教程

    热门讨论 2009-08-04 18:33:50
    运用地统计学进行空间分析基本包括以下几个步骤,即数据探索性分析,空间连续性的量化模型,未知点属性值的估计,对未知点局部及空间整体不确定性的预测。
  • 半方差函数与Moran_sI在土壤微量元素空间分布研究中的应用_以寿光市为例.pdf
  • 把离散半方差模型投资组合问题,推广到连续时间情形.引进恰当的状态约束,将原问题转化为一个有约束的随机最优控制问题.利用经典Lagrange理论,将其进一步转化为无约束随机LQ(Li near Quadratic)最优控制问题.进而借助...
  • 基于加权半方差的含风电电力系统旋转备用效益研究.pdf
  • 已经实现二维样点的空间结构变异性分析编码,可分析样点的同向及异向半方差、协方差,绘制方差图,也可作为空间插值的数据预处理
  • 主要介绍了基于python计算滚动方差(标准差)talib和pd.rolling函数差异详解,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 14,245
精华内容 5,698
关键字:

半方差