CN110287601B - 一种毛竹胸径年龄二元联合分布精确估算方法 - Google Patents

Abstract

本发明公开了一种毛竹胸径年龄二元联合分布精确估算方法，包括以下步骤：一，将毛竹林胸径年龄按株数汇总成二维统计数据；二，绘制毛竹胸径年龄的三维频数直方图；三，得若干种常用二元Copula密度函数的参数向量；四，绘制毛竹胸径年龄的若干种二元Copula密度函数图；五，选取预备二元Copula密度函数图；六，得毛竹胸径年龄最优二元Copula密度函数模型；七，得毛竹胸径年龄二元Copula分布估算函数，进而实现毛竹胸径年龄二元联合分布的精确估算。本发明具有根据胸径与年龄的边缘分布值即可建立毛竹胸径年龄的二元Copula分布(密度)函数；根据相邻两期毛竹连续清查数据，可推算出中间年份胸径年龄的联合分布(密度)值的优点。

Description

一种毛竹胸径年龄二元联合分布精确估算方法

技术领域

本发明涉及一种区域森林生物量精确估算模型的建立方法，特别是一种毛竹胸径年龄二元联合分布精确估算方法。

背景技术

毛竹是我国南方重要的一种森林类型，也是最重要、最典型的竹林资源类型，中国是竹林资源最丰富的国家之一。毛竹也是一种特殊的地域性异龄林植物，其胸径生长与高生长在第一年内完成，且毛竹林多有隔年出笋成竹即大小年习性，故毛竹林林分年龄结构按2年划分1个龄级，因此同时描述毛竹林的胸径与年龄能更完整地反映其资源的结构特征。据研究，毛竹林分还具有很强的碳汇功能(周国模，姜培坤，2004；方晰等，2002；李意德等，1998；阮宏华等，1997)，是未来重要的生物质能源材料，因此如何精确估算区域尺度毛竹生物量显得非常重要。

测树因子的二元概率分布对区域森林生物量精确估算及现代森林经营活动起着非常重要的作用，它可以精确预估树木各径阶、各龄阶株数与蓄积变化、生长量测定、产量预估、枯损量计算等，还对抚育间伐的科学设计与森林碳库精确估算提供科学依据。到目前为止，二元广义β分布和Johnson's SBB分布被广泛应用于研究测树因子的胸径、树高二元联合分布(金星姬，2013；Li et al.，2002；Wang and Rennolls，2007；Hafley et al.，1976；H.T.SCHREUDER et al.，1977；FASHENG LI,2002；Mingliang Wang，2007)，二元Weibull分布与二元Beta分布函数描述树冠内细枝叶面积分布的二维结构(Temesgenetal.,2003)，二元Weibull分布与二元最大熵函数用来描述毛竹胸径、年龄的联合分布(葛宏立，2008；刘恩斌，2010)。

已有研究测树因子二元联合分布的方法都有如下缺陷：1)边缘分布模型必须相同：二维Beta分布、二维Weibull分布、二维Sbb的边缘分布必须是Beta分布、Weibull分布与Sbb分布，但在实际应用时，当两个测树因子的边缘分布不同时，常用二元分布函数不能描述这两个测树因子的联合分布，因此，常用二元分布模型的适用性不广；2)不是由边缘分布确定对应的联合分布：常用研究测树因子二元联合分布的方法是先建立其联合分布函数，后再推导联合分布对应的边缘分布(刘恩斌，2010；葛宏立，2008)，故建立的联合分布函数不是由对应边缘分布函数确定的，事实上边缘分布函数的类型与理论较联合分布函数更多、更完善，但理论与应用价值远没有联合分布函数大，因此如何根据测树因子边缘分布建立其二元联合分布的方法有非常重要的意义；3)不能根据相邻两期毛竹连续清查样地推算中间年份的数据：毛竹连续清查样地每5年进行一次复查，常用测树因子二元联合分布的参数拟合需给定毛竹胸径与年龄的联合分布数据，而只根据相邻两期连续清查数据是不能确定中间任一年份的毛竹胸径年龄联合分布数据，故常用二元联合分布函数不能精确推算出中间4年的毛竹连续清查数据，因此不能准确预估省域尺度毛竹生物量的动态变化。

随着科学研究的不断推进，在数学研究领域快速发展的Copula理论(Nelsen,2006)为解决上述相关非正态变量的建模问题提供了新的思路。Copula理论的核心思想就是任何一个二元联合分布函数都可拆分为一个可包含相关性大小和相关结构信息的Copula函数和两个对应的边缘分布函数，因此Copula理论在建模时的灵活性和适用性方面有着无可比拟的优势，根据常用研究测树因子二元联合分布模型的诸多不足，将Copula理论应用到林业与生态这一领域是当下林业研究中的一大趋势，如何应用Copula函数建立毛竹胸径、年龄二元联合分布的模型更是目前林业研究中亟需解决的一大难题。

发明内容

本发明的目的在于，提供一种基于Copula函数建立毛竹胸径年龄二维分布模型的方法。它具有1)、根据胸径与年龄的边缘分布值即可建立毛竹胸径年龄的二元Copula分布(密度)函数；2)、根据相邻两期毛竹连续清查数据，可推算出中间年份胸径年龄的联合分布(密度)值的优点。

本发明的技术方案：一种毛竹胸径年龄二元联合分布精确估算方法，包括以下步骤：

步骤一，将研究区域N个毛竹样地的连续清查数据按株数汇总成毛竹林胸径年龄二维统计数据；

步骤二，根据步骤一的二维统计数据分别计算毛竹胸径的边缘分布值u与年龄边缘分布值v，以u与v为自变量，用matlab软件绘制毛竹胸径年龄的三维频数直方图；

步骤三，以毛竹胸径与年龄累计分布值为自变量，采用极大似然法，用若干种常用二元Copula密度函数分别拟合步骤一的二维统计数据，得各二元Copula密度函数的参数向量；

步骤四，用matlab软件分别绘制若干种毛竹胸径年龄二元Copula密度函数图；

步骤五，将步骤四中所有毛竹胸径年龄二元Copula密度函数图与步骤二中得到的三维频数直方图进行对比，从中选取与三维频数直方图形状接近的几幅Copula密度函数图作为预备二元Copula密度函数图；

步骤六，根据AIC信息准则与确定系数公式，分别计算若干种毛竹胸径年龄二元Copula密度函数的AIC值与确定系数R²，选取AIC值最小、R²最大且与步骤二中的三维频数直方图形状最接近的二元Copula密度函数模型作为毛竹胸径年龄最优二元Copula密度函数模型D(u,v)；

步骤七，根据概率论中的密度函数与分布函数的转换关系，可得毛竹胸径年龄二元Copula分布估算函数，进而实现毛竹胸径年龄二维分布的精确估算。

前述的一种毛竹胸径年龄二元联合分布精确估算方法中，所述步骤六中AIC的值计算式为：

式中：D(u_i,v_i；θ)为所选二元Copula密度函数；θ为所选二元Copula函数参数向量；m为所选二元Copula密度函数的参数个数；L为所选二元Copula密度函数的极大似然函数；N为要拟合的原始数据的样本量，当所选模型为二维Gaussian Copula、Frank Copula、Gumbel Copula、Clayton Copula函数时，m等于1，所选模型为二维t Copula函数时，m等于2，u_i,v_i(i＝1，2，…N)分别为原始样本数据统计得到的边缘分布值。

前述的一种毛竹胸径年龄二元联合分布精确估算方法中，所述步骤六中确定系数R²的计算式为：

式中y_i为因变量的第i个实测值，

为因变量的第i个估计值，N为样本个数，

为因变量的平均值。

前述的一种毛竹胸径年龄二元联合分布精确估算方法中，所述步骤三、步骤四与步骤六中若干种二元Copula密度函数包括二维正态Copula密度函数、二维t-Copula密度函数、二维Gumbel Copula密度函数、二维Clayton Copula密度函数和二维Frank Copula密度函数。

前述的一种毛竹胸径年龄二元联合分布精确估算方法中，所述二元Copula密度函数图的绘制具体包括以下步骤：

a、通过步骤一中的统计数据计算毛竹胸径与年龄的边缘分布值；

b、以毛竹胸径与年龄的边缘分布值为自变量，采用极大似然法，分别拟合5种二元Copula密度函数参数向量，进而由毛竹胸径与年龄的边缘分布值估算出毛竹胸径年龄二元Copula密度函数值；

c、以毛竹胸径与年龄的边缘分布值为自变量，毛竹胸径年龄二元Copula密度函数值为因变量，用matlab软件绘制毛竹胸径年龄的二元Copula密度函数图。

与现有技术相比，本发明针对常用二元分布(密度)函数的不足，采用二元Copula函数描述毛竹胸径年龄的二元分布，基于Copula函数建立毛竹胸径年龄二维分布模型的方法，适合于胸径与年龄的任意边缘分布，具体应用时，不需要确定边缘分布函数的类型，根据实测数据估算出胸径与年龄的边缘分布函数值即可；具体地，毛竹胸径年龄二元Copula函数的建立首先毛竹样地的连续清查数据确定胸径与年龄的边缘分布值，采用AIC信息准则计算毛竹胸径与年龄的边缘分布值在各个预备二元Copula密度函数模型下的AIC值，最终选取AIC值最小的函数模型作为最优Copula函数模型，拟合最优Copula函数模型的参数α，得到最终的二元Copula概率密度估算函数。

二元Copula函数的参数拟合只需胸径与年龄边缘分布实测值，不需胸径年龄实测数据的联合分布(密度)值，根据相邻两期毛竹连续清查数据，即可推算出中间年份胸径年龄的联合分布(密度)值；二元Copula分布(密度)函数可以描述测树因子间的非线性相关性，进而刻画出测树因子不同部位的相关性。

综上，本发明具有1)、根据胸径与年龄的边缘分布值即可建立毛竹胸径年龄的二元Copula分布(密度)函数；2)、根据相邻两期毛竹连续清查数据，可推算出中间年份胸径年龄的联合分布(密度)值的优点。

附图说明

图1是毛竹胸径年龄累计分布三维频数直方图；

图2是二维正态Copula联合密度函数图；

图3是二维t-Copula联合密度函数图；

图4是二维Gumbel Copula联合密度函数图；

图5是二维Clayton Copula联合密度函数图；

图6是二维Frank Copula联合密度函数图；

图7是毛竹胸径年龄二元Copula分布函数图；

图8是毛竹胸径年龄二维概率密度函数图。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明，但并不作为对本发明限制的依据。

实施例：一种毛竹胸径年龄二元联合分布精确估算方法，包括以下步骤：

步骤一，将研究区域N个毛竹样地的连续清查数据按株数汇总成毛竹林胸径年龄二维统计数据，浙江省于1979年建立了森林资源连续清查体系，以5年为一个复查周期，共设置固定样地4250个，样点格网为4km×6km，样地形状为正方形，边长28.28m，面积800m²。利用2009年毛竹177个连续清查样地进行研究，毛竹样地的基本情况为：每个样地毛竹株数22～897不等，样地内毛竹胸径为5～15cm，年龄1～4度以上(当年生竹记为1度竹；2-3年生竹记为2度竹，依此类推)。样地的调查因子包括土层厚度、坡向(1：北，2：东北，3：东，4：东南，5：南，6：西南，7：西，8：西北，9：无坡向即山顶)、坡位(1：脊部，2：上部，3：中部，4：下部，5：谷部，6：平地)、坡度、海拔(10m-1200m)、样地平均胸径、样地株数及样地内每株毛竹的胸径与年龄等。本实施例中样地的选取是全省区域的宏观模型，所以将245个样地的数据合并统计到一个表上，见表1。表中密度值为某径阶、某龄阶的毛竹株数，分布值为从起始径阶、起始龄阶到某径阶、某龄阶的毛竹株数累积值；

表1浙江省毛竹林胸径、年龄样地统计数据

步骤二，根据步骤一的二维统计数据分别计算毛竹胸径的边缘分布值u与年龄边缘分布值v，以u与v为自变量，用matlab软件绘制毛竹胸径年龄累计分布的三维频数直方图，如图1所示；

步骤三，以毛竹胸径与年龄累计分布值为自变量，采用极大似然法，用若干种常用二元Copula密度函数分别拟合步骤一的二维统计数据，得各二元Copula密度函数的参数向量；

步骤四，用matlab软件分别绘制若干种毛竹胸径年龄二元Copula密度函数图，常用Copula函数有正态Copula函数、t-Copula函数、阿基米德Copula函数，其中GumbelCopula函数、ClaytonCopula函数、Frank Copula函数为最常用的3种阿基米德Copula函数(张蕾，2017)，因此，本发明中采用这5种二元Copula函数对毛竹胸径年龄的联合密度进行估计，得到的二元Copula密度函数图分别如图2-图6所示；

步骤五，将步骤四中所有毛竹胸径年龄二元Copula密度函数图与步骤二中得到的三维频数直方图进行对比，从中选取与三维频数直方图形状接近的几幅Copula密度函数图作为预备二元Copula密度函数图，从图1与图2-图6的对比中可看出，二维Gumbel Copula与二维t Copula联合密度图与图1类似，其它二元Copula函数与图1有较大差异，因此优先选择Gumbel Copula与t Copula密度函数图作为预备二元Copula密度函数图；

步骤六，根据AIC信息准则与确定系数公式，分别计算若干种毛竹胸径年龄二元Copula密度函数的AIC值与确定系数R²，选取AIC值最小、R²最大且与步骤二中的三维频数直方图形状最接近的二元Copula密度函数模型作为毛竹胸径年龄最优二元Copula密度函数模型D(u,v)，根据AIC的值计算式，计算二维t Copula、二维Gumbel Copula概率密度函数的AIC值与确定系数R²分别为-14.3104、-19.5196与0.9676、0.9841，所以二维Gumbel Copula概率密度函数描述毛竹胸径年龄二元联合分布最好，为最优Copula函数模型；

对所有二元Copula密度函数图用AIC信息准则再次进行验算，Gaussian Copula的AIC值为2，t Copula的AIC值为-14.3104，Clayton Copula的AIC值为-1.2248，FrankCopula的AIC值为-2.4059，Gumbel Copula的AIC值为-19.5196；确定二维Gumbel Copula概率密度函数为最优Copula函数模型；

步骤七，根据概率论中的密度函数与分布函数的转换关系，可得毛竹胸径年龄二元Copula分布估算函数，进而实现毛竹胸径年龄二维分布的精确估算，此时毛竹胸径年龄联合密度函数的估计值，见下表2：

表2省域尺度毛竹实测概率及相应二元正态Copula概率密度函数的估计值

从表2可以看出，二维Gumbel Copula概率密度函数对胸径年龄联合密度的估计精度很高，在下尾部的估计值相对实际值略高。

所述步骤六中AIC的值计算式为：

式中：D(u_i,v_i；θ)为所选二元Copula密度函数；θ为所选二元Copula函数参数向量；m为所选二元Copula密度函数的参数个数；L为所选二元Copula密度函数的极大似然函数；N为要拟合的原始数据的样本量，当所选模型为二维Gaussian Copula、Frank Copula、Gumbel Copula、Clayton Copula函数时，m等于1，所选模型为二维t Copula函数时，m等于2，u_i,v_i(i＝1，2，…N)分别为原始样本数据统计得到的边缘分布值。

所述步骤六中确定系数R²的计算式为：

式中y_i为因变量的第i个实测值，

为因变量的第i个估计值，N为样本个数，

为因变量的平均值。

所述步骤三、步骤四与步骤六中若干种二元Copula密度函数包括二维正态Copula密度函数、二维t-Copula密度函数、二维Gumbel Copula密度函数、二维Clayton Copula密度函数和二维Frank Copula密度函数。

所述二元Copula密度函数图的绘制具体包括以下步骤：

a、通过步骤一中的统计数据计算毛竹胸径与年龄的边缘分布值；

b、以毛竹胸径与年龄的边缘分布值为自变量，采用极大似然法，分别拟合5种二元Copula密度函数参数向量，进而由毛竹胸径与年龄的边缘分布值估算出毛竹胸径年龄二元Copula密度函数值；

具体通过(1)式，由毛竹胸径与年龄的边缘分布函数值算出毛竹胸径、年龄的二元Copula分布函数值，并由matlab软件绘制相应的图，如图7所示；再根据(2)式，由毛竹胸径、年龄的二元Copula分布函数值推导出相应的二元Copula密度函数值；最后根据(3)式，由二元Copula密度函数值可得毛竹胸径、年龄二维密度函数值；

c、以毛竹胸径与年龄的边缘分布值为自变量，毛竹胸径年龄二元Copula密度函数值为因变量，用matlab软件绘制毛竹胸径年龄的二元Copula密度函数图，如图8所示；

所述(1)式为H(x,y)＝C(F(x),G(y))，其中H(x,y)为二维随机变量(x，y)的联合分布，其边缘分布分别为F(x),G(y)；

所述(2)式为Copula函数C(u,v)的概率密度函数D(u,v)，

u＝F(x),v＝G(y)；

所述(3)式为H(x,y)的概率密度函数h(x,y)，

为变量x的边缘概率密度函数，

为变量y的边缘概率密度函数；

其中，(2)式与(3)式的关系为：

式中D(u,v)为二元Copula密度函数值，h(x,y)为毛竹胸径、年龄二维密度函数值。

Copula函数最早由Sklar(1959)提出：将多个一维分布连接成多维分布的函数，同时每个边缘分布必须是均匀分布。

Copula函数基本理论：根据定义(Nelsen,2006)，二元Copula函数定义如下：

(1)C(·,·)在每个维度空间的定义域均为[0,1]；

(2)C(·,·)有零基面且二维递增的；

(3)对于任意的变量u,v∈[0,1],，均满足C(u,1)＝u,C(1,v)＝v；

同时，对于定义域内的任意一点(u,v)，都有0≤C(u,v)≤1。

根据定义可知，Copula函数有如下性质：

1、边界条件：对任意二元变量(u,v)∈[0,1]×[0,1]，均满足C(u,0)＝0,C(0,v)＝0，C(u,1)＝u,C(1,v)＝v，即只要有一个变量为0，Copula函数值就为0，只要有一个变量为1，Copula函数值由另一个变量决定。

2、递增性：Copula函数在定义域内严格单调非递减，对于定义域内的变量u₁,v₁,u₂,v₂，且u₁≤u₂,v₁≤v₂，则C(u₂,v₂)-C(u₂,v₁)-C(u₁,v₂)+C(u₁,v₁)≥0，即如果u,v∈[0,1]的值同时增大，则Copula函数值一定是非减少的，同理如果一个自变量不变，则Copula函数值随着另一个变量值的增大而增大或不变。

3、Frechet边界：任一个Copula分布都有Frechet上界C⁺(u,v)＝max(u+v-1,0)与下界C-(u,v)＝min(u,v)，即C-(u,v)≤C(u,v)≤C⁺(u,v)。

Sklar定理：令二维随机变量(x，y)的联合分布为H(x,y)，边缘分布分别为F(x),G(y)，则存在着一个Copula函数C，(1)式：

H(x,y)＝C(F(x),G(y))

若F(x)和G(y)是连续的，则C是唯一的，反过来，对任何概率分布F(x)和G(y)和Copula函数C，则式(1)定义的H(x,y)一定是一个联合分布函数，其边缘为F(x)和G(y)。

设D(u,v)是Copula函数C(u,v)的概率密度函数，由(1)式得(2)式：

令u＝F(x),v＝G(y)，则

为变量x的边缘概率密度函数，

为变量y的边缘概率密度函数，则H(x,y)的概率密度函数h(x,y)为(3)式：

(3)式表明，一个联合概率密度函数h(x,y)拆解成两部分，前一部分D(x,y)为Copula密度函数，能精确完整地描述随机变量之间的相关结构，后一部分f(x)g(y)则为边缘概率密度函数的乘积。

毛竹胸径年龄的联合分布对毛竹林的经营与区域尺度生物量的精确估算意义重大，而Copula函数种类较多，适应性广，在林业与生态中有非常广泛的应用，结合本文的研究，得出如下结论：

1)基于Copula函数建立毛竹胸径年龄二维分布模型的方法，适合于胸径与年龄的任意边缘分布，具体应用时，不需要确定边缘分布函数的类型，根据实测数据估算出胸径与年龄的边缘分布函数值即是Copula函数的自变量。

2)毛竹胸径年龄二元Copula函数的建立分为两大步：第一步确定胸径与年龄的边缘分布值，第二步确定最优Copula函数及拟合相应Copula函数的参数。

3)二维Gumbel Copula分布(密度)函数对毛竹胸径年龄联合分布的测量精度最高，R²＝0.9841，为最优Copula函数。

二元Copula分布(密度)函数与毛竹胸径年龄常用二元概率分布(密度)函数的相比较有着以下不同：

一、两者自变量不同：毛竹胸径年龄常用二元概率分布(密度)函数，如二元Weibull、二元Sbb、二元Beta函数，其自变量为胸径与年龄，因变量为胸径年龄的联合分布(密度)值，而二元Copula分布(密度)函数的自变量为毛竹胸径与年龄的边缘分布函数值，因变量为Copula连接函数值。

二、拟合分布(密度)函数参数所需条件不同：常用二元概率分布(密度)函数参数的拟合必需胸径年龄实测数据的联合分布(密度)值，而二元Copula分布(密度)函数的参数拟合只需要毛竹胸径与年龄实测数据的边缘分布值，因此拟合常用二元概率分布(密度)函数参数所需的条件比二元Copula分布(密度)函数高的多，这也从侧面说明二元Copula分布(密度)函数的应用更广、实用性更强。

三、胸径与年龄边缘分布函数的类型：应用二元Beta分布(密度)函数拟合胸径年龄联合分布时，必须要求胸径与年龄的边缘分布函数类型为一元Beta分布(密度)函数，二元Weibull、二元Sbb函数也有类似的要求，如果胸径与年龄边缘分布函数的类型不同，则常用二元分布(密度)函数不适用，而二元Copula分布(密度)函数适用于任意类型的胸径与年龄边缘分布。

四、联合分布函数的建立：通常由测树因子的联合分布是可以确定对应的边缘分布，边缘分布是很难确定对应联合分布的，因此常用二元联合分布(密度)函数不是由对应的边缘分布(密度)函数确定的，不管在理论还是实践中，测树因子的边缘分布的类型众多且较容易确定，但联合分布的价值较边缘分布更大，因此如何根据测树因子的边缘分布确定对应的联合分布意义重大，以上的研究表明，二元Copula分布(密度)函数可以把联合分布与对应的边缘分布连接起来，从而估算出毛竹胸径年龄的二元联合分布(密度)值，进而克服常用二元联合分布(密度)函数的不足。

五、模型的检验：实测胸径年龄的联合分布(密度)值是拟合常用二元概率分布(密度)函数参数必须要的数据，而我们用二元分布(密度)函数的最终目的是为了得到胸径、年龄的联合分布(密度)的估计值，故为了验证模型的可靠性与稳定性，实测数据需分为建模样本与检验样本，以对常用二元概率分布(密度)函数进行参数的拟合与检验，而二元Copula分布(密度)函数参数的拟合只需要胸径与年龄的边缘分布函数值(不需要实测胸径、年龄的联合分布(密度)值)，因此不需要检验样本对模型进行检验。

六、区域尺度毛竹胸径年龄联合分布(密度)预估：浙江省于1979年开始进行森林资源连续清查，且每5年进行一次，到目前为止共完成进行了7次资源的连续清查，如果要对全省毛竹资源进行时间上的预测，则7期的数据远远不够，这时必须要按年份对相邻两期数据进行加密，而用常用二元概率分布(密度)函数进行预测时，首先要拟合函数参数，而函数参数的拟合必须要有毛竹胸径年龄的联合分布(密度)值，而根据相邻两期毛竹连续清查数据，不能推算出中间年份胸径年龄的联合分布(密度)值，即常用二元概率分布(密度)函数不能对森林连续清查数据进行加密，而毛竹胸径与年龄的生长很特殊，有别于其他植物，根据相邻两期胸径与年龄实测数据的边缘分布(密度)值是可以推算出中间年份胸径年龄的边缘分布(密度)值，故二元Copula分布(密度)函数可以对相邻两期毛竹连续清查数据进行加密，从而可精准预测区域尺度毛竹胸径年龄联合分布(密度)，进行对区域尺度的毛竹生物量进行预估。

各二元Copula分布(密度)函数的比较：

对比图2、图3、图6，二维正态Copula、二维t Copula与二维Frank Copula密度函数具有对称的尾部，它们无法描述毛竹胸径与年龄非对称的尾部相关性，其中二维t Copula函数有较厚的尾部，对测树因子尾部相关性较为敏感，能更好地描述测树因子的对称性尾部相关性，二元Clayton Copula与二元Gumbel Copula函数具有不对称的尾部，能描述毛竹胸径与年龄非对称的尾部相关性，其中二元Gumbel Copula函数的图像呈J字形，表明该函数对毛竹胸径与年龄的上尾部变化较为敏感，能更好地描述测树因子上尾高、下尾低的尾部相关性，二元Clayton Copula函数的密度函数图像呈L形，表明该函数对测树因子下尾部的变化较敏感，能更好地描述测树因子上尾低、下尾高的尾部相关性。

Claims (5)

1.一种毛竹胸径年龄二元联合分布精确估算方法，其特征在于，包括以下步骤：

步骤一，将研究区域N个毛竹样地的连续清查数据按株数汇总成毛竹林胸径年龄二维统计数据；

步骤二，根据步骤一的二维统计数据分别计算毛竹胸径的边缘分布值

与年龄边缘分布值

，以

与

为自变量，用matlab软件绘制毛竹胸径年龄的三维频数直方图；

步骤三，以毛竹胸径与年龄累计分布值为自变量，采用极大似然法，用若干种常用二元Copula密度函数分别拟合步骤一的二维统计数据，得各二元Copula密度函数的参数向量；

步骤四，用matlab软件分别绘制若干种毛竹胸径年龄二元Copula密度函数图；

步骤五，将步骤四中所有毛竹胸径年龄二元Copula密度函数图与步骤二中得到的三维频数直方图进行对比，从中选取与三维频数直方图形状接近的几幅Copula密度函数图作为预备二元Copula密度函数图；

步骤六，根据AIC信息准则与确定系数公式，分别计算若干种毛竹胸径年龄二元Copula密度函数的AIC值与确定系数R²，选取AIC值最小、R²最大且与步骤二中的三维频数直方图形状最接近的二元Copula密度函数模型作为毛竹胸径年龄最优二元Copula密度函数模型

；

步骤七，根据概率论中的密度函数与分布函数的转换关系，可得毛竹胸径年龄二元Copula分布估算函数，进而实现毛竹胸径年龄二元分布的精确估算。

2.根据权利要求1所述的一种毛竹胸径年龄二元联合分布精确估算方法，其特征在于：所述步骤六中AIC的值计算式为：

式中：

为所选二元Copula密度函数；

为所选二元Copula函数参数向量；m为所选二元Copula密度函数的参数个数；L为所选二元Copula密度函数的极大似然函数；N为要拟合的原始数据的样本量，当所选模型为二元 Gaussian Copula、Frank Copula、Gumbel Copula、Clayton Copula函数时，m 等于1，所选模型为二元t Copula函数时，m 等于2，

，分别为原始样本数据统计得到的边缘分布值。

3.根据权利要求1所述的一种毛竹胸径年龄二元联合分布精确估算方法，其特征在于：所述步骤六中确定系数R²的计算式为：

式中

为因变量的第i个实测值，

为因变量的第i个估计值，N为样本个数，

为因变量的平均值。

4.根据权利要求1所述的一种毛竹胸径年龄二元联合分布精确估算方法，其特征在于：所述步骤三、步骤四与步骤六中若干种二元Copula密度函数包括二元正态Copula密度函数、二元t-Copula密度函数、二元Gumbel Copula密度函数、二元Clayton Copula密度函数和二元Frank Copula密度函数。

5.根据权利要求4所述的一种毛竹胸径年龄二元联合分布精确估算方法，其特征在于：所述二元Copula密度函数图的绘制具体包括以下步骤：

a、通过步骤一中的统计数据计算毛竹胸径与年龄的边缘分布值；

b、以毛竹胸径与年龄的边缘分布值为自变量，采用极大似然法，分别拟合5种二元Copula密度函数参数向量，进而由毛竹胸径与年龄的边缘分布值估算出毛竹胸径年龄二元Copula密度函数值；

c、以毛竹胸径与年龄的边缘分布值为自变量，毛竹胸径年龄二元Copula密度函数值为因变量，用matlab软件绘制毛竹胸径年龄的二元Copula密度函数图。

Images