回归估计

2024-10-10

回归估计（精选9篇）

回归估计篇1

考虑线性回归模型

$Y = X β + ε, ε ~ (0, σ^{2} Ι_{n}) (1)$

式 (1) 中Y是n×1观测变量, X是n×p的列满秩矩阵 (即R (X) =p) , β是p×1的未知参数, ε是n×1随机误差。线性回归模型中未知参数β的最小二乘估计 (LS) 一直以来倍受重视。但当设计的矩阵病态时, LS估计的精度很差, 为了克服这种弊端, β的有偏估计受到了统计学家的关注, 例如岭估计。在文献[1]中, K.J.LIU提出了一种新的有偏估计——LIU估计:

$\hat{β} (d) = (X^{'} X + Ι)^{- 1} (X^{'} X + d Ι) \hat{β} 。$

式中 $\hat{β}$ 是β的LS估计。

研究LIU估计的文献有很多, 其中文献[2]研究了广义LIU估计及其优良性, 指出了LIU估计是线性模型最小二乘估计的有偏估计, 它是可容许估计且局部优于最小二乘估计。文献[3,4,5]分别讨论了不同模型下各种估计的影响分析问题。本文中我们主要研究LIU估计的影响分析问题, 得到不同模型下LIU估计的影响度量式。

本文中考虑下面的模型:

数据删除模型:

$Y_{[i]} = X^{'}_{[i]} β + ε_{[i]}, ε_{[i]} ~ (0, σ^{2} Ι_{n - 1}) (2)$

方差扰动模型:

Y=Xβ+ε, ε～ (0, σ2G-1) (3)

均值漂移模型:

$Y = X β + D γ + ε, ε ~ (0, σ^{2} Ι_{n}) (4)$

为了度量在不同模型对LIU估计的影响程度, 类似cook距离, 我们定义广义cook距离

$D (Μ, c) = (\hat{β}_{d} (f) - \hat{β}_{d})^{'} Μ (\hat{β}_{d} (f) - \hat{β}_{d}) / c (5)$

式 (5) 中, $\hat{β}$ d (f) 代表模型式 (3) , 式 (4) , 式 (5) 的LIU估计, $\hat{β}$ d代表原模型式 (1) 的LIU估计。

引理1[6]:设

$A = (\begin{array}{l} A_{11} \\ A_{21} \end{array} \begin{array}{l} A_{12} \\ A_{22} \end{array})$

可逆。若 $| A_{11} | \neq 0$ , 则

$\begin{array}{l} A^{- 1} = (\begin{array}{l} A_{11} \\ A_{21} \end{array} \begin{array}{l} A_{12} \\ A_{22} \end{array})^{- 1} = \\ (\begin{array}{l} A_{11}^{- 1} + A_{11}^{- 1} A_{12} A_{22.1}^{- 1} A_{21} A_{11}^{- 1} \\ - A_{22.1}^{- 1} A_{21} A_{11}^{- 1} \end{array} \begin{array}{l} - A_{11}^{- 1} A_{12} A_{22.1}^{- 1} \\ A_{22.1}^{- 1} \end{array}) 。 \end{array}$

若 $| A_{22} | \neq 0$ , 则

$\begin{array}{l} A^{- 1} = (\begin{array}{l} A_{11} \\ A_{21} \end{array} \begin{array}{l} A_{12} \\ A_{22} \end{array})^{- 1} = \\ (\begin{array}{l} A_{11.2}^{- 1} \\ - A_{22}^{- 1} A_{21} A_{11.2}^{- 1} \end{array} \begin{array}{l} - A_{11.2}^{- 1} A_{12} A_{22}^{- 1} \\ A_{22}^{- 1} + A_{22}^{- 1} A_{21} A_{11.2}^{- 1} A_{12} A_{22}^{- 1} \end{array}) 。 \end{array}$

其中:A22.1 =A22 -A21 A $_{11}^{- 1}$ A12 , A11.2 =A11 -A12 ×A $_{22}^{- 1}$ A21 。

引理2[6]:若A与 (A+BCD) 均可逆, 则

(A+BCD) -1=A-1-A-1B (I+CA-1B) CA-1。

定理1:模型式 (2) 的LIU估计 $\hat{β}$ d (i) 与 $\hat{β}$ d有如下关系:

$\hat{β}_{d} (i) = \hat{β}_{d} - \frac{(X^{'} X + Ι)^{- 1} (X^{'} X + d Ι) (X^{'} X)^{- 1} x_{i} {\hat{e}}^{'}_{i}}{1 - h_{i i}} +$

$\begin{array}{l} \frac{(X^{'} X + Ι)^{- 1} x_{i} x^{'}_{i} \hat{β}_{d}}{1 - \bar{h}_{i i}} + \frac{(X^{'} X + Ι)^{- 1} x_{i} (h_{i i} y^{'}_{i} - x^{'}_{i} \hat{β})}{(1 - \bar{h}_{i i}) (1 - h_{i i})} - \\ \frac{(X^{'} X + Ι)^{- 1} x_{i} x^{'}_{i} (X^{'} X + Ι)^{- 1} (X^{'} X + d Ι) (X^{'} X)^{- 1} x_{i} {\hat{e}}^{'}_{i}}{(1 - \bar{h}_{i i}) (1 - h_{i i})} \\ (6) \end{array}$

$\begin{array}{l} 式 (6) 中 \bar{h}_{i i} = x^{'}_{i} (X^{'} X + Ι)^{- 1} x_{i}, h_{i i} = x^{'}_{i} (X^{'} X)^{- 1} x_{i}, \\ {\hat{e}}^{'}_{i} = y^{'}_{i} - x^{'}_{i} \hat{β} 。 \end{array}$

且若在式 (6) 中取M= (X′X+I) , c=σ2, 则有

$\begin{array}{l} D (Μ, c) = \frac{(\hat{β}_{d} (i) - \hat{β}_{d})^{'} Μ (\hat{β}_{d} (i) - \hat{β}_{d})}{c} = \\ [- \frac{\hat{e}_{i} x^{'}_{i} (X^{'} X)^{- 1} (X^{'} X + d Ι)}{1 - h_{i i}} + \frac{{\hat{β}}^{'}_{d} x_{i} x^{'}_{i}}{1 - \bar{h}_{i i}} + \\ \frac{(h_{i i} y_{i} - {\hat{β}}^{'} x_{i}) x^{'}_{i}}{(1 - \bar{h}_{i i}) (1 - h_{i i})} - \\ {\hat{e}_{i} x^{'}_{i} (X^{'} X)^{- 1} (X^{'} X + d Ι) (X^{'} X + Ι)^{- 1} x^{'}_{i} x_{i}} \times \\ {(1 - \bar{h}_{i i}) (1 - h_{i i})}^{- 1}] \times \\ \frac{(X^{'} X + Ι)^{- 1}}{σ^{2}} [- \frac{(X^{'} X + d Ι) (X^{'} X)^{- 1} x_{i} {\hat{e}}^{'}_{i}}{1 - h_{i i}} + \\ \frac{x_{i} x^{'}_{i} \hat{β}_{d}}{1 - \bar{h}_{i i}} + \frac{x_{i} (h_{i i} y^{'}_{i} - x^{'}_{i} \hat{β})}{(1 - \bar{h}_{i i}) (1 - h_{i i})} - \\ \frac{x_{i} x^{'}_{i} (X^{'} X + Ι)^{- 1} (X^{'} X + d Ι) (X^{'} X)^{- 1} x_{i} {\hat{e}}^{'}_{i}}{(1 - \bar{h}_{i i}) (1 - h_{i i})}] (7) \end{array}$

证明:根据LIU估计的定义可得到模型式 (2) 的LIU估计

$\begin{array}{l} \hat{β}_{d} (i) = (X^{'}_{[i]} X_{[i]} + Ι)^{- 1} (X^{'}_{[i]} X_{[i]} + \\ d Ι) (X^{'}_{[i]} X_{[i]})^{- 1} X^{'}_{[i]} Y_{[i]} (8) \end{array}$

注意到:

X′[i]X[i]=X′X-xix′i, X′[i]Y[i]=X′Y-xiy′i,

(X′[i]X[i]+dI) = (X′X+dI) -xix′i (9)

根据引理2, 可以得到

(X′[i]X[i]+I) -1= (X′X+I-xix′i) -1= (X′X+I) -1+ (X′X+I) -1xi (1-x′i (X′X+

I) -1xi) -1x′i (X′X+I) -1 (10)

(X′[i]X[i]) -1= (X′X-xix′i) -1= (X′X) -1+

(X′X) -1xi (1-x′i (X′X) -1xi) -1×x′i (X′X) -1 (11)

将式 (9) 、式 (10) 、式 (11) 代入式 (8) 并整理, 可以得出式 (6) , 在式 (5) 中将 $\hat{β}$ d (f) 用 $\hat{β}$ d (i) 代替并取M= (X′X+I) , c=σ2, 可以推出式 (7) , 证毕。

定理2:模型式 (3) 的LIU估计 $\hat{β}$ d (G) 与 $\hat{β}$ d有如下关系

$\begin{array}{l} \hat{β}_{d} (G) = \hat{β}_{d} + (X^{'} X + Ι)^{- 1} (X^{'} X + d Ι) (X^{'} X)^{- 1} X^{'} \bar{G} (Ι - Η \bar{G})^{- 1} \hat{e} + (X^{'} X + Ι)^{- 1} X^{'} \bar{G} (Ι - Η_{1} \bar{G})^{- 1} X (\hat{β}_{d} - \hat{β}) + (X^{'} X + Ι)^{- 1} X^{'} \bar{G} (Ι - Η_{1} \bar{G})^{- 1} (Η - Η_{2}) \bar{G} (Ι - \\ Η \bar{G})^{- 1} \hat{e} (12) \end{array}$

$\begin{array}{l} 式 (12) 中 ∶ \hat{β}_{d} = (X^{'} X + Ι)^{- 1} (X^{'} X + d Ι) (X^{'} X)^{- 1} X^{'} Y ‚ \\ \hat{β} = (X^{'} X)^{- 1} X^{'} Y ‚ Η = X (X^{'} X)^{- 1} X^{'} ‚ Η_{1} = X (X^{'} X + Ι)^{- 1} \times X^{'} ‚ Η_{2} = X (X^{'} X + Ι)^{- 1} (X^{'} X + d Ι) (X^{'} X)^{- 1} X^{'} ‚ \hat{e} = Y - X \hat{β} ‚ \bar{G} = Ι - G 。 \end{array}$

若在式 (5) 中取M= (X′X+I) , c=σ2, 则有

$\begin{array}{l} D (Μ, c) = \frac{(\hat{β}_{d} (G) - \hat{β}_{d})^{'} Μ (\hat{β}_{d} (G) - \hat{β}_{d})}{c} = \\ [{\hat{e}}^{'} (Ι - Η \bar{G})^{- 1} \bar{G} X (X^{'} X)^{- 1} (X^{'} X + d Ι) + (\hat{β}_{d} - \\ \hat{β})^{'} X^{'} (Ι - Η_{1} \bar{G})^{- 1} \bar{G} X + {\hat{e}}^{'} (Ι - Η \bar{G})^{- 1} \bar{G} (Η - \\ Η_{2}) (Ι - Η \bar{G})^{- 1} \bar{G} X] \frac{(X^{'} X + Ι)^{- 1}}{σ^{2}} [(X^{'} X + \\ d Ι) (X^{'} X)^{- 1} X^{'} \bar{G} (Ι - Η \bar{G})^{- 1} \hat{e} + X^{'} \bar{G} (Ι - \\ Η_{1} \bar{G})^{- 1} (Η - Η_{2}) \bar{G} (Ι - Η \bar{G})^{- 1} \hat{e} + \\ X^{'} \bar{G} (Ι - Η_{1} \bar{G})^{- 1} X (\hat{β}_{d} - \hat{β})] (13) \end{array}$

证明:由LIU估计的定义得到模型式 (3) 的LIU估计

$\begin{array}{l} \hat{β}_{d} (G) = (X^{'} G X + Ι)^{- 1} (X^{'} G X + \\ d Ι) (X^{'} G X)^{- 1} X^{'} G Y (14) \end{array}$

记 $\bar{G} = Ι - G ‚ Η = X (X^{'} X)^{- 1} X^{'} ‚ Η_{1} = X (X^{'} X + Ι)^{- 1} X^{'}$ , 由引理2得到

$\begin{array}{l} (X^{'} X + Ι - X^{'} \bar{G} X)^{- 1} = (X^{'} X + Ι)^{- 1} + (X^{'} X + Ι)^{- 1} \\ X^{'} \bar{G} (Ι - X (X^{'} X + Ι)^{- 1} X^{'} \bar{G})^{- 1} X (X^{'} X + Ι)^{- 1} X (15) \\ (X^{'} X - X^{'} \bar{G} X)^{- 1} = (X^{'} X)^{- 1} + (X^{'} X)^{- 1} X^{'} \bar{G} (Ι - \\ X (X^{'} X)^{- 1} X^{'} \bar{G})^{- 1} X (X^{'} X)^{- 1} (16) \end{array}$

注意到

$\begin{array}{l} Ι + (Ι - Η \bar{G})^{- 1} Η \bar{G} = (Ι - Η \bar{G})^{- 1}, \\ Ι + (Ι - Η_{1} \bar{G})^{- 1} Η_{1} \bar{G} = (Ι - Η_{1} \bar{G})^{- 1} (17) \end{array}$

将式 (15) 、式 (16) 、式 (17) 代入式 (14) 并整理得到式 (12) , 在式 (5) 中将 $\hat{β}$ d (f) 用 $\hat{β}$ d (G) 代替得到式 (13) , 证毕。

定理3:模型式 (4) 的LIU估计 $\hat{β}$ d (γ) 与 $\hat{β}$ d有如下关系

$\begin{array}{l} \hat{β}_{d} (γ) = \hat{β}_{d} + d (X^{'} X + Ι_{p})^{- 1} (X^{'} X)^{- 1} X^{'} D (Ι_{r} - \\ D^{'} X (X^{'} X)^{- 1} X^{'} D)^{- 1} D^{'} \hat{e} + (X^{'} X + Ι_{p})^{- 1} X^{'} D (2 Ι_{r} - \\ D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D)^{- 1} D^{'} X (\hat{β}_{d} - \hat{β}) - (d - 1) \times \\ (X^{'} X + Ι_{p})^{- 1} X^{'} D (2 Ι_{r} - D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D)^{- 1} \times D^{'} X (X^{'} X + Ι_{p})^{- 1} (X^{'} X)^{- 1} X^{'} D (Ι_{r} - D^{'} X (X^{'} X)^{- 1} \times X^{'} D)^{- 1} D^{'} \hat{e} + (X^{'} X + Ι_{p})^{- 1} X^{'} D (2 Ι_{r} - D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D)^{- 1} [D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D - (d + 1) Ι_{r}] \times (Ι_{r} - D^{'} X (X^{'} X)^{- 1} X^{'} D)^{- 1} D^{'} \hat{e} (18) \end{array}$

若在式 (5) 中取M= (X′X+Ip) , c=σ2, 则有

$\begin{array}{l} D (Μ, c) = \frac{(\hat{β}_{d} (γ) - \hat{β}_{d})^{'} Μ (\hat{β}_{d} (γ) - \hat{β}_{d})}{c} = \\ {d {\hat{e}}^{'} D (Ι_{r} - D^{'} X (X^{'} X)^{- 1} X^{'} D)^{- 1} D^{'} X (X^{'} \times X)^{- 1} + (\hat{β}_{d} - \hat{β})^{'} X^{'} D (2 Ι_{r} - D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D)^{- 1} D^{'} X - (d - 1) {\hat{e}}^{'} D (Ι_{r} - \\ D^{'} X (X^{'} X)^{- 1} X^{'} D)^{- 1} D^{'} X (X^{'} X)^{- 1} (X^{'} X + Ι_{p})^{- 1} X^{'} D (2 Ι_{r} - D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D)^{- 1} D^{'} X + {\hat{e}}^{'} D (Ι_{r} - D^{'} X (X^{'} X)^{- 1} X^{'} D)^{- 1} [D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D - (d + 1) Ι_{r}] (2 Ι_{r} - D^{'} X (X^{'} X + Ι_{p})^{- 1} D^{'} X)^{- 1} X^{'} D} \frac{(X^{'} X + Ι_{p})^{- 1}}{σ^{2}} {d (X^{'} \times X)^{- 1} X^{'} D (Ι_{r} - D^{'} X (X^{'} X)^{- 1} \times X^{'} D)^{- 1} D^{'} \hat{e} + X^{'} D (2 Ι_{r} - D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D)^{- 1} \times \\ D^{'} X (\hat{β}_{d} - \hat{β}) - (d - 1) X^{'} D (2 Ι_{r} - D^{'} X (X^{'} X + Ι_{p})^{- 1} \times X^{'} D)^{- 1} D^{'} X (X^{'} X + Ι_{p})^{- 1} (X^{'} X)^{- 1} X^{'} D (Ι_{r} - D^{'} X (X^{'} X)^{- 1} X^{'} D)^{- 1} D^{'} \hat{e} + X^{'} D (2 Ι_{r} - D^{'} X (X^{'} X + \\ Ι_{p})^{- 1} X^{'} D)^{- 1} [D^{'} X (X^{'} X + Ι_{p})^{- 1} X^{'} D - (d + 1) Ι_{r}] \times (Ι_{r} - D^{'} X (X^{'} X)^{- 1} X^{'} D)^{- 1} D^{'} \hat{e}} (19) \end{array}$

证明:记

$C = (X, D), α = (\begin{array}{l} β \\ γ \end{array})$

, 模型式 (4) 可写为

$Y = C α + ε (20)$

由LIU估计的定义得到模型式 (20) 的LIU估计为

$\hat{β}_{d} (γ) = (C^{'} C + Ι)^{- 1} (C^{'} C + d Ι) (C^{'} C)^{- 1} C^{'} Y (21)$

将C= (X, D) 代入C并将其写成矩阵形式得到

$\begin{array}{l} \hat{β}_{d} (γ) = (\begin{array}{l} X^{'} X + Ι_{p} \\ D^{'} X \end{array} \begin{array}{l} X^{'} D \\ D^{'} D + Ι_{r} \end{array})^{- 1} \times \\ (\begin{array}{l} X^{'} X + d Ι_{p} \\ D^{'} X \end{array} \begin{array}{l} X^{'} D \\ D^{'} D + d Ι_{r} \end{array}) \times \\ (\begin{array}{l} X^{'} X \\ D^{'} X \end{array} \begin{array}{l} X^{'} D \\ D^{'} D \end{array})^{- 1} (\begin{array}{l} X^{'} \\ D^{'} \end{array}) Y (22) \end{array}$

应用引理1将 $\hat{β}$ d (γ) 展开并整理得到式 (18) , 在式 (5) 中将 $\hat{β}$ d (f) 用 $\hat{β}$ d (γ) 代替得到式 (19) , 证毕。

参考文献

[1]Liu KJ.Anewclass of biased estimate in linear regression.CommIn Stat, 1993;22 (3) :393—402

[2]陈德英, 张尚立.广义LIU估计及其优良性, 科学技术与工程, 2008;8 (12) :3272—3274

[3]林路.协方差阵扰动模型岭估计的影响分析.工程数学学报, 1995, 12 (3) :83—88

[4]王铭, 石磊.协方差分析模型的影响分析.云南大学学报, 2003;25 (5) :391—394

[5]张尚立, 覃红.约束条件下线性模型协方差阵扰动的影响分析.数学物理学报, 2006;26A (4) :621—628

[6]王松桂.线性模型理论及其应用.合肥:安徽教育出版社, 1987

回归估计篇2

应用支持向量回归估计预测陀螺误差系数

针对目前小样本容量陀螺误差系数预测精度不高的.问题,本文将支持向量回归估计引入到陀螺误差系数的预测研究中.通过对某型陀螺某项误差系数的预测,并且对比分析该方法与目前通用的AR模型预测方法的预测效果,结果表明本文采用的支持向量回归估计具有更高的预测精度.

作者：焦巍王宏力刘光斌 JIAO Wei WANG Hong-li LIU Guang-bin 作者单位：第二炮兵工程学院,西安,710025 刊名：电光与控制 ISTIC PKU英文刊名：ELECTRONICS OPTICS & CONTROL 年，卷(期)： 13(5) 分类号：V241.5 关键词：支持向量机支持向量回归估计误差系数预测 AR模型

回归估计篇3

1 贝叶斯估计和传统估计的比较

贝叶斯法与传统方法的区别主要在于传统方法做出的推断基于频率的基础，对估计值的统计性质实际上是对大量重复实验做出的推断，而不能对当前估计值的统计性质给出一个明确的答案，在实际应用中含有更多的不确定性。

经典统计学认为参数是某个未知的常数，可以通过样本信息来获得。样本是来自于某个总体，而该总体用分布函数或密度函数来描述，总体分布往往含有未知参数，如位置参数、形状参数等。用样本信息构造的统计量来估计总体参数的过程就称为参数估计。

而贝叶斯统计学认为参数的取值是随机的，经验的或数据的信息形成参数的先验概率分布。在得到样本之前，从先验分布中产生了某个样本，即参数值。在样本已知的情况下，似然函数包含了所有的样本信息，先验分布和似然函数的乘积得到参数的后验分布。该参数的估计使得损失函数的数学期望最小。可以看出，在无先验信息或弱先验信息的情况下，贝叶斯估计与极大似然估计理论上应一致，在可靠的先验信息情况下，贝叶斯估计的误差应该更小。

2 置信区间与贝叶斯区间估计的区别

经典统计学用样本信息构造的统计量来估计总体参数，为了衡量估计值的误差范围，需要给出区间估计。区间估计往往通过统计量的分布来获得。上下限通常是某个统计量。由于统计量是随机变量，由抽样差异导致每次计算的置信区间都是不同的，这些区间要么包含总体参数，要么不包含，也就是说区间是随机的。理论上，当抽样次数达到很大时，大约有100×（1-a）个区间覆盖真值。1-a为置信水平，置信水平越高，区间越窄。

贝叶斯区间估计则大大不同，由于假设参数是随机的，前述当获得参数的后验分布后，可以计算出该参数落入任何一个区域内的概率，可以给定概率，比如0.95，通常可以找到一些区域，贝叶斯区间选择使得该区间上的后验密度值不少于相同概率下其他区域上的值。也就是最高后验密度的区间即相同概率下最短的区间。由此得到的区间可以解释为该参数落入其中的概率为1-a。这与实际中的理解吻合，而避免了重复抽样的牵强的解释。

3 模拟实验结果分析

该模拟实验通过对多因变量的回归模型的模拟结果，清楚地揭示了先验信息在贝叶斯分析中的作用，以及贝叶斯回归分析及传统回归分析的区别和联系（见表1）。模拟数据来自模型：

在弱先验信息下（先验估计值的变异程度为33.3），给定两组不同的初值，iterations=1000，从以上表1可以看出，在弱先验信息的条件下，MCMC得到的结果与传统最小二乘法得到的结果基本一致；两者在beta项上的估计都不准确，但是在自变量参数值的估计上比较准确，从SD可以看出最小二乘法得到的估计要更精确。当调整先验信息的把握度，（先验估计值的变异程度为0.033)MCMC结果与真实的结果一致，比传统方法更为准确。

由此可以得出当先验信息准确且有把握的情况下，调整先验信息在MCMC的结果中所占权重，结果偏好于先验信息。MCMC得到的估计其实是在经验与数据之间的折衷。同时，可以通过对旧的数据得到的参数估计作为新的先验信息，更新得到新的模型。

4 CEAC曲线解释的比较

该部分以弱先验信息下绘制的成本-效果接受曲线为例来说明两种方法在参数解释上的差异。引入意愿支付法后，自然得到：

INB(Rc)=Rc*增量效果-增量成本=Rc*△E-△C

Rc表示增量效果意愿支付的货币值，INB代表增量净效益，prob代表INB>0的概率。

必须在得到的联合分布后，从而得到增量净效益的分布，得到增量净效益>0的概率。用联合分布可以反应出△E、△C之间的联系，当两者独立时，联合分布等于各自分布的乘积。

图1中Rc代表患者对增加单位疗效而愿意支付的钱，代表选择治疗方案二和治疗方案一相比所获得的收益。图中在点Rc=3.9时，p=0.5，贝叶斯派的解释为当增加单位效果患者愿意付超过3.9元时，则选择治疗有较大的把握能够获得盈利，而不是亏本。这里的把握程度会随着Rc的取值慢慢增大，即概率值增大。但是按频率派的解释，当患者愿意支付超过3.9元时，重复100次这样的治疗，大约有超过一半的治疗会盈利，而不是亏本，至于该次治疗的增量-净效益会不会>0，是未知的，无法做出推断。

同时从图1中可以看出，强先验信息下，对本次数据结果的INB>0的概率要大于弱先验信息下的概率，表明在有可靠先验信息的情况下，对同一数据的结果所得INB>0的把握度（信念）是不一样的。在信息准确的情况下，△E、△C分布越集中，同一Rc下INB>0的把握越大。进一步说明，先验信息（经验的、数据的）将影响最终的决策，甚至可以得出相反的结论。例如，当Rc=11.36时，弱先验下的概率为0.55，强先验下的概率为0.877。弱先验下就有可能选择第二种方案，而强先验下只可能选择第一种方案。

5 结果与讨论

本文主要探讨了贝叶斯方法和频率法之间的区别和联系，并通过模拟实验给出了一个清晰的结果。在此基础上分析了采用贝叶斯方法在成本—效果分析中的作用，计算了成本—效果接受曲线，比较了不同先验信息下的差异以及两种方法在解释上的差异，阐明了在国外经济学分析中偏好贝叶斯方法的原因。当然本文对成本—效果分析中的一些问题讨论的还不够全面，诸如，成本—效果可接受曲线的一般计算过程、方法等。另外还可以将研究扩展到对增量成本效果比的区间估计做比较，同时还可以在bootstrap法与贝叶斯方法得到的区间估计中做一个作为方法间比较的另一证据。

参考文献

[1]张玉哲,吴晶,孙利华.药物经济学评价中的贝叶斯法介绍[J].中国药物经济学,2009(2):51-56.

[2]阿诺德泽尔纳,张尧庭,蒋传海,等.阿诺德泽尔纳[M].上海:上海财经大学出版社,2005:7.

[3]O'Hagan,A.and J.W.Stevens.The probability of costeffectiveness[J].BMC Med Res Methodol,2002(2):5.

[4]朱慧明,林静.贝叶斯计量经济模型[M].北京:科学出版社,2009:66-67.

[5]叶慈南,曹伟丽.应用数理统计[M].上海:机械工业出版社,2007:42-67.

回归估计篇4

本文考虑纵向数据下半参数回归模型:yij=xij′β+g(tij)+eij,i=1,…,m,j=1,…,ni.基于最小二乘法和一般的非参数权函数方法给出了模型中参数β,回归函数g(・)和误差方差σ2的.估计,并在适当条件下证明了估计量的r(r≥2)阶平均相合性.

作者：田萍薛留根 TIAN Ping XUE Liu-gen 作者单位：田萍,TIAN Ping(许昌学院数学系,河南,许昌,461000)

薛留根,XUE Liu-gen(北京工业大学应用数理学院,北京,100022)

回归估计篇5

我国是一个农业大国, ‘民以食为天’, 农业发展在我国国民经济中占有特别重要的地位。随着农业经济的发展加快, 农民收入不断增加, 农民消费水平不断提高, 农村市场的巨大潜力在不断释放, 都最终会影响到整个国民经济的各个方面。因此, 我们必须把发展农业始终放在国民经济发展的首位, 进一步加强农业的基础地位, 走可持续发展的道路, 增加农民收入, 提高农民消费能力, 从粗放型转变为集约型, 发展循环经济, 为国民经济的发展提供更加长久的动力。

二、指标描述

经济增长指一个国家或地区生产的物质产品和服务的增加, 它关心的重点是物质方面的进步、生活水平的提高, 衡量经济增长最重要的指标为GDP。农业增加值反映了农业生产经营活动的最终成果和对社会的贡献, 是农业发展的最主要指标, 用AGR表示。本文对各地区生产总值和农业增加值2006年末的横截面数据进行研究, 获得数据如下表。

三、经典回归分析

以GDP纵坐标, AGR为横坐标, 画出二者的散点图如图1。

从图中可以看出, 地区生产总值和农业增加值间有正相关关系, 且广东、山东和江苏三省的数据明显高于其它地区的数据, 究其原因, 这三省处于沿海地区, 经济比较发达, 农业现代化较好, 而且由于其地理位置优越, 因而农业发展较好, GDP也较高。这种数据的存在必将影响经典的最小二乘回归分析的结果, 但是这三省都是经济大省, 故分析时不能随意将其删除, 需要对其进行基于异常值的回归分析。下面我们对其分别进行经典的最小二乘回归分析和稳健回归分析, 并对二者进行比较。

用Qstat软件分析, 进行经典回归分析结果拟合方程如上表

由软件输出结果F=34.323, 对应的p值为0.0000, 模型显著, 但调整决定系数为0.5262, 即自变量AGR只能解释因变量GDP的52.62%, 模型拟合效果不佳。

由Qstat软件输出的诊断图如图2

由学生化残差图可以看出存在异常值, 致使拟合的回归直线向上偏离大部分数据。

四、稳健回归

从上面的分析中可以看出, 当数据中存在异常值时, LS估计的结果显著变坏。但是, 基于异常值的稳健M估计受异常值的影响不是非常明显, 即起到抗异常值的作用, 下面进行稳健回归.

由QSTAT输出的结果可拟合回归方程为:

由F=254.47, 回归模型显著。调整后的决定系数为0.8942, 即自变量AGR可以解释因变量GDP的89.42%, 模型拟合效果良好。

由Qstat软件输出的稳健回归的诊断图如图3。

有诊断图可以看出, 基于稳健估计拟合的模型并不受异常值的影响而使回归曲线向上偏移, 而是将异常值的影响减小。

五、结论及建议

1. 由上面的比较可以看出, 本文采用经典回归农业增加值只能解释GDP的52.62%, 若采用M估计农业增加值则能解释GDP的89.42%, 因而M估计性能更好一些。

2. 农业发展对经济增长的贡献作用巨大。从稳健估计得出的回归模型可以看出, 农业增加值每增加1个单位, GDP增加5.62个单位。虽然随着经济的发展第二产业和第三产业的比重会扩大, 但农业的发展仍然是GDP增长的主要拉动因子, 继续大力发展农业, 是保证国民经济稳定发展的重要举措。

3. 中国目前正处于高速发展的阶段, 通过农业发展对经济增长的贡献的分析, 亦为如何更快促进中国经济发展带来一定的启示。即制度创新、技术创新使农业经济增长的方式从粗放型转向集约型, 加大农村人力资本的投入, 努力提高农村劳动力的素质, 农业循环经济发展之路, 实现农业可持续发展, 使农业在经济增长中充分发挥其不可替代的基础作用。

摘要：虽然农业在经济增长中的贡献份额日渐减少, 但是农业发展对经济增长的作用却并未减少, 本文通过对31省、直辖市经济增长与农业发展的数据进行经典回归和基于异常值的稳健回归分析比较, 研究农业发展对经济增长的贡献。

关键词：经济增长,农业发展,稳健,回归

参考文献

[1]王斌会:经济管理模型的多变量统计方法及分析系统TSTAT[M].中国统计出版社, 2005

[2]林卿等:可持续农业经济发展论—经济可持续发展论丛[M].中国环境科学出版社, 2002

[3]乔玲:浅谈新形势下农业发展对策[J].农业经济, 2006年第5期

回归估计篇6

统计过程控制 (statistics process control, SPC) 是监控制造过程异常, 辅助维护过程稳定进而保证产品质量的有效手段。然而, 传统SPC方法单纯采用控制图发现产品质量特性的异常波动, 而不能诊断出关于异常的具体信息[1], 如常用的均值控制中, 检测到均值发生异常时, 偏移幅度究竟有多大无法直接获知, 而这种信息对异常诊断和过程改进都极其重要。针对这一问题, 学术界开展了将机器学习方法与传统SPC相结合进行过程异常智能诊断的研究。目前, 研究最多的是基于人工神经网络 (artificial neural network, ANN) 的方法[2-4]。

近年来, 基于统计学习理论的支持向量机 (support vector machine, SVM) [5]被提出并得到迅速发展, 其在学习能力、泛化能力以及实用性方面比ANN更具有优越性。在SPC领域也有了SVM的应用研究, 其中, 支持向量回归 (supportvector regression, SVR) 作为SVM在回归领域的延伸, 开始被用于对过程异常幅度的预测和估计。如Cheng等[6]提出用径向基核函数型SVR估计过程均值偏移量, 并通过仿真实验证实了该方法比ANN方法估计精度更高。

按照SVR原理, 原始空间非线性分布的数据通过核函数的作用被映射到高维特征空间中, 变得近似服从线性分布, 故而可以进行线性回归[7]。因此, 核函数是决定SVR性能的关键因素, 且不同的应用情况本质上要求有与之相匹配的核函数, 才能实现最合理的空间映射。然而, 在实际应用中, 核函数的构造和选取至今没有完善的理论依据, 常用的几种基本核函数具有不同映射特性, 针对不同应用场合表现出较大性能差异。只采用单个核对未知分布的样本进行处理并不合理, 因此, 近年来出现了多核学习方法的研究[8-10]。多核的基本思想即是将多种基本核函数进行集成, 以发挥各自的优势, 以增强决策函数的可解释性[11]。

为更加准确地对制造过程中监控的质量特性均值变量的偏移幅度进行估计, 以提高异常在线诊断能力, 本文提出一种采用粒子群优化 (PSO) 算法对多个简单核函数进行优化组合而实现的多核SVR (PSOMK-SVR) , 它通过PSO算法对各核参数、组合权重系数以及SVR的惩罚系数进行联合优化而得到, 再用于对过程均值偏移的回归估计。

1 支持向量回归数学模型

SVR是对SVM基本原理的推广, 用于解决非线性回归估计问题。SVR解决的问题可描述为[12]:对样本集T= { (xi, yi) |i=1, 2, …, l}, 其中, xi∈Rn, 为n维输入变量, yi∈R, 为xi对应输出值, l为样本数, 寻求样本最优拟合函数f ∈F= {f|f:Rn→ R}, 使之对非样本输入也能精确估计输出。针对该问题, SVR通过非线性映射φ:Rn→ Rm (m ≥n) 将原始的样本数据从输入空间 (n维) 映射到高维特征空间 (m维) , 使样本数据转化为易于进行线性拟合的高维数据, 再构造出线性拟合超平面:

式中, ω 为超平面法向量;b为超平面偏置量;〈·, ·〉表示特征空间中的内积运算。

SVR寻求最优拟合超平面等效于求解SVM的最大间隔问题, 并引入松弛变量使拟合具有一定的柔性, 得到如下优化问题:

其中, ξi和ξi*是表征系统上下限的松弛变量;ε为不敏感系数;C是取值为正常数的惩罚因子, 用于在拟合精度与泛化性能间取得平衡;|y- 〈ω, φ (x) 〉-b|是ε线性不敏感损失函数, 定义如下:

引入拉格朗日乘子αi、αi*、λi、λi*, 利用Wolfe对偶技巧, 建立原问题的拉格朗日对偶问题:

在式 (4) 中引入核函数K (xi, xj) = 〈xi·xj〉, 并利用KKT条件算出参数αi、αi*、b值, 最后得到SVR回归拟合函数为

2 多核学习原理及PSOMK-SVR构建方法

2.1 核函数类型及性质

从SVR的基本原理可知, 核函数是决定其回归性能的关键因素。核函数是满足Mercer条件的对称函数, 几种常用基本核函数如表1所示。

表1中, u和v为输入空间向量。根据核函数的定义可推出其具有如下基本性质[12]:① 若K1和K2是两个核, α1和α2是两个正实数, 则K =α1K1+α2K2也是核;② 若K1和K2是两个核, 则K = K1K2也是核;③ 若K1是核, 则K1的指数K (u, v) =exp (K1) 也是核。

2.2 多核学习原理及凸组合核

从本质上讲, 核函数对内积具有什么样的计算 (运算形式以及系数取值) , 就会形成什么样的特征空间。选择一种核函数即确定了从样本空间到特征空间的映射形式, 也就决定了SVR不同的回归拟合精度以及泛化能力。一般说来, 线性核具有全局特性, 相距很远的数据点都能影响到核函数值, 易获得较好泛化能力;径向基核具有局部特性, 只有相距很近的数据点对核函数值才有影响, 具有较好的学习能力。而一个学习模型性能的好坏通常由学习能力和泛化能力两方面决定, 因此可考虑对多种核函数进行集成, 综合其优势。前述核函数基本性质为此提供了思路, 多核学习正是利用这些性质集成多个核函数, 其中, 最直接的多核构建方法即将基本核函数按下式进行凸组合:

其中, K*为组合核函数, Ki、ωi分别为第i核函数及其权重, 满足ωi∈[0, k], 。

由此构建的多核函数综合了各基本核的特征映射能力, 将不同核函数形成的多个特征空间进行组合, 使数据在新的特征空间中得到更好表达[11]。

2.3 利用粒子群算法构建优化多核函数

本文按照式 (6) 构建多核函数, 选取表1所示3种基本核进行组合, 现要解决的首要问题就是确定组合权重系数。此外, 还需确定各基本核函数自身核参数的取值。核参数很大程度上决定了基本核发挥其映射特性的程度, 例如径向基核的局部映射能力与其核宽度参数γ具有密切联系, 不同参数可能导致映射结果出现极大差异。同时, SVR的惩罚系数C也是决定其泛化性能的关键参数。要获得性能良好的多核SVR, 需要这几类参数同时取得合适的值。为此, 本文引入粒子群优化算法[13]对这几类参数进行联合优化, 构建多核SVR。

PSO算法具有参数少、效率高等优点, 它以粒子Xi代表d维空间里具有记忆能力的点Xi= (xi1, xi2, …, xid) 。迭代过程中每个粒子根据之前的种群最优和个体最优位置信息调整自身位置, 调整公式为

式中, r1、r2为区间[0, 1]上的随机数;Vi(k)、Xi(k)分别为Xi第k次迭代时的速度和位置;Pi(k)为Xi第k次迭代中的最优位置;Pg(k)为k次迭代中所有粒子的最优位置;ω为惯性权重因子;c1、c2为加速度因子。

本文的PSO粒子编码为 (λ, coef, γ, C, ω1, ω2) 。λ∈ [0, 1.0], coef ∈ [0, 3.0], γ∈ [0, 1], 分别为多项式核和径向基核参数 (为简化计算, 多项式核固定为3次多项式, 即d=3) ;C∈ [1, 200];ω1∈ [0, 3]、ω2∈ [0, 3], 分别为线性核和多项式核权重, 径向基核权重ω3=3-ω1-ω2。为准确评价多核SVR的泛化估计能力, 引入五折交叉验证, 并采用决定系数R2作为精度指标。R2为估计值与真实值相关系数的平方, 在[0, 1]之间取值, 越接近1表明估计精度越高, 定义如下[13]:

式中, l为测试样本个数;yi为样本i真实值;^yi为样本i估计值。

利用PSO构建优化多核函数的过程如图1所示。

3 基于PSOMK-SVR的过程均值异常监测模型

为将PSOMK-SVR用于制造过程均值异常幅值估计, 将累积和 (CUSUM) 控制图与

该监测模型由两大模块构成:①用于训练出PSOMK-SVR结合, 构建如图2所示的监测模型。其中, CUSUM图基于序贯分析原理, 具有累积效应, 对过程均值的微小偏移表现得比图更灵敏, 利于更早发现过程异常, 也更早估计出偏移幅值。

PSOMK-SVR的离线学习模块;②制造过程均值变量异常在线监测模块。离线学习包括如下步骤:

(1) 根据采集的大量过程变量历史数据统计出其分布参数 (本文考察的均值变量, 根据中心极限定理其近似服从正态分布, 分布参数为均值μ和标准差σ) 。

(2) 以蒙特卡洛仿真方法生成训练与测试样本。幅值估计的原始样本向量定义为Xs= (xt*-w/2, xt*-w/2+1, …, xt*-1, xt*, xt*+1, …, xt*+w/2-2, xt*+w/2-1) T, 其中w为观测窗口宽度, t*代表CUSUM控制图检出偏移异常的观测点, Xs包括以t*为中心的连续w个观测值 (t*位于窗口中央) 。对应每个样本向量, 以发生偏移异常时变量幅值的实际偏移量作为目标值。

(3) 对原始样本向量依次进行标准化、特征融合与归一化处理得到训练与测试样本。标准化处理使得训练出的学习器独立于具体过程参数, 具有普遍适用性;特征融合能更好地刻画样本特征, 提高估计精度, 为此提取了经标准化处理的原始样本向量的6种统计特征[14] (均值 (mean) 、标准差 (SD) 、偏度 (skewness) 、峰度 (kurtosis) 、均方值 (mean-square value) 、自相关量 (autocorrela-tion) ) , 并与标准化处理的原始样本向量串联融合;归一化处理将数据统一到同一尺度范围内, 利于提高计算效率。最终计算样本向量为:X′t*= (x′t*-w/2, x′t*-w/2+1, …, xt*-1, x′t*, xt*+1, …, x′t+*w/2-2, x′t+*w/2-1, x′mean, …, x′autocorr) T。

(4) 随机选取部分训练样本, 利用前述PSO方法对参数 (包括核参数、权重系数和SVR惩罚系数) 进行优化, 最后以全部训练样本训练出有效的PSOMK-SVR, 并以测试样本集进行测试。

在线异常监测包括两个环节, 首先由CUSUM控制图对均值进行监控, 即通过在线提取新的均值计算CUSUM统计量, 与控制限值比对以判断是否出现异常。CUSUM统计量计算如下[15]:

其中, Ci+和Ci-为上下单侧累积和, 初值C0+=C0-=0;K为参考值;μ0为均值期望。当Ci+或Ci-超过控制限H时, 判定过程异常。然后由PSOMK-SVR构成的偏移幅值估计环节在过程出现失控后, 再等待半个窗口宽度的观测点, 提取出原始样本向量, 并经过与生成训练样本时相同的标准化、特征融合及归一化处理后, 输入到PSOMK-SVR中估计出当前过程均值的偏移量值。

4 仿真实验验证

仿真实验在MATLAB7.11 软件中进行, SVR采用LIBSVM工具箱。硬件环境为CPU2.2GHz, 内存2GB, 操作系统为WindowsXP3.2。

4.1 实验数据生成

利用蒙特卡洛仿真生成训练与测试样本, 仿真模型如下:

X (t) =μ0+r (t) σ0+bδσ0 (11)

其中, X (t) 为t观测点的均值变量值, μ0为均值变量期望, σ0为均值变量标准差, r (t) ~ N (0, 1) ;δ表示均值偏移幅值为σ0倍数, b为开关变量, 在异常前后分别取值0和1。

设μ0=80, σ0=5, δ∈{-3.5, -3.0, …, -0.5, 0.5, …, 3.0, 3.5}, 将δ的每种取值代入式 (11) 中生成长度为1000的数据流, 且设定在t=25时发生异常, 即使b从0变为1。CUSUM控制图参数定义如表2所示, 将生成的数据流输入到CUSUM控制图中, 当检测到异常时, 以当前数据点xt*为中心提取观测窗口w =8个均值数据构成原始样本向量。然后对原始样本向量进行标准化、特征融合及归一化处理得到训练和测试样本向量。生成训练样本共10×14=140个, 测试样本共10×14=140个。

4.2 PSOMK-SVR生成

SVR不敏感系数取为经验值, ε=0.01;PSO粒子取值范围为:λ∈ [0, 1.0], coef ∈ [0, 3.0], γ∈ [0, 1], C ∈ [1, 200], ω1∈ [0, 3], ω2∈ [0, 3];PSO的种群规模为20, 最大迭代次数为100, 加速度因子c1=c2=1.5;为更好地平衡全局搜索和局部搜索能力, 惯性权重按如下二次曲线动态调整[13]:

式中, ωstart=0.9, 为初始值;ωend=0.4, 为结束值;T为当前迭代次数;Tmax为最大迭代次数。

迭代得最优参数λ=0.72, coef =0.74, γ=1.0, C=126, ω1=1.86, ω2=0.075, 迭代过程如图3所示。代入获得的最优参数, 并以全部训练样本作为输入训练出PSOMK-SVR。

4.3 结果与讨论

为验证PSOMK-SVR估计过程均值偏移幅度的有效性和准确性, 将其与BP神经网络及通常认为效果最好的径向基核函数 (RBF) 型SVR (RBF-SVR) 进行对比。其中, BP网络为单隐层结构, 隐层激发函数tansig, 输出层激发函数logsig, 最大迭代次数1000, 训练精度0.001, 学习速率0.1, 采用有动量加自适应学习速率的梯度下降法训练;RBF-SVR通过网格搜索和五折交叉验证法进行参数整定。三者都用生成的训练样本集进行训练, 并作用于相同的测试样本集。比较的指标包括对全部测试样本估计的均方误差ΔMSE、决定系数R2, 以及每种偏移幅度的平均绝对百分误差ΔMAPE, 结果如图4所示。

图4a~图4c分别给出了三者对全部测试样本的估计情况, 从图上直观看出, PSOMK-SVR的估计值与真实值的重合度最高, 根据指标ΔMSE和R2量值看来, PSOMK-SVR的ΔMSE值最小, R2值最大。图4d中给出了ΔMAPE的对比情况, 可清晰地看出PSOMK-SVR的估计误差明显比BP神经网络的估计误差小, 也小于RBF-SVR, 尤其在偏移幅度本身较小情况下优势更明显。可见, SVR在小样本条件下回归估计性能明显好于神经网络;经过对多种核函数进行参数及权重系数寻优得到的优化多核SVR相对于单个核函数的SVR, 回归估计性能得到进一步提高。

5 应用验证

为进一步验证PSOMK-SVR的有效性和实用性, 以它为基础开发了过程异常在线监测软件, 对某型汽车变速箱输出一挡齿轮齿坯内孔 (图5) 精车工序进行应用验证。该工序属齿轮自动加工线的组成部分, 由车削中心以内圆车刀加工完成, 并通过线边检测台进行在线测量。根据前期累积的大量历史数据统计得到该齿坯孔径尺寸均值, 标准差, 以此作为过程基本参数建立图和CUSUM控制图。鉴于PSOMK-SVR的生成过程已经对训练样本作了标准化处理, 因此可将训练好的模型直接集成到监测软件中使用。

软件以容量为5的采样子组对该尺寸测量值进行获取, 并计算出子组均值记录到图中, 如图6所示, 同时通过CUSUM图进行异常判别。程序在第52号观测点即被检测出向下偏移异常,

在等待4个观测点 (即半个观测窗口宽度) 后, 提取49~56号观测点数据到PSOMK-SVR中, 进行偏移幅度估计, 得到δ=-1.74, 与子组均值标准差相乘得到估计的偏移量 Δ =-0.0012mm。为检验软件的估计效果, 没有立即对过程进行干预, 继续加工直到获取更多的X值 (一共获取了62组观测值) 。然后计算52号观测点到62号观测点间共11个X值的均值, 并与值作差, 得到实际发生偏移量 Δ′ =-0.0013mm。Δ与 Δ′ 相对误差为8.3%, 由此证明PSOMK-SVR能对实际制造过程变量的偏移幅值进行较为准确的估计, 具有一定实用性。

6 结论

传统SPC控制图不能直接提供过程变量发生异常的幅度, 不利于异常的诊断及过程的改进。通过机器学习方法对过程异常幅度进行智能估计能有效解决这一问题, 其中, 基于ANN以及基于单核函数SVR的过程异常参数估值方法已被证明是有效手段。

本文提出一种对多个不同特性的基本核函数进行凸组合构建SVR, 并通过PSO对包括核函数参数、组合权重系数和SVR惩罚系数的全部参数进行联合优化的方法。该参数寻优过程以五折交叉验证计算训练样本的决定系数均值作为PSO的粒子适应度值, 将寻优目标导向为使SVR具有更好的拟合精度和泛化能力。

回归估计篇7

在轻型舰载鱼雷发射的过程控制和参数分析中, 动力学模型往往比较复杂, 如何利用有限的实验数据, 进行鱼雷发射过程中一些主要参数的分析及快速预报是十分有意义的。在潜艇鱼雷发射装置管内运动参数设计中, 首先要弄清鱼雷管峰压的影响因素及其函数关系。文献[1]利用基于最小二乘法的回归分析方法, 建立了发射管内峰压Pm的一元线性回归方程, 并进行了显著性检验。

偏最小二乘回归是一种新的多元数据分析方法, 集多元线性回归、典型相关分析和主成分分析的基本功能为一体, 将建模预测类型的数据分析方法与非模型式的数据认识方法有机地结合起来, 使得模型精度、稳健性、实用性得到提高[2,3,4], 已广泛应用于各领域的研究中[5,6,7,8,9,10]。本文应用偏最小二乘回归方法, 进一步探索鱼雷水面发射时影响鱼雷管峰压Pm的主要因素, 建立回归方程并进行误差分析。

1偏最小二乘回归模型

1.1建模原理

设有p自变量X={x1, x2, …, xp}和q个因变量Y={y1, y2, …, yp}, 观测n个样本点, 分别在X与Y中提取出成分t1和u1 (即t1是x1, ..., xp的线性组合, u1是y1, ..., yq的线性组合) , 要求t1和u1应尽可能大地携带他们各自数据表中的变异信息, 同时t1与u1的相关程度能够达到最大使得t1和u1应尽可能好的代表数据表X和Y, 同时自变量的成分t1对因变量的成分u1又有最强的解释能力。

在第一个成分t1和u1被提取后, 分别实施X对t1的回归以及Y对u1的回归。如果回归方程已经达到满意的精度, 则算法终止;否则, 将利用X被t1解释后的残余信息以及Y被t1解释后的残余信息进行第二轮的成分提取。如此往复, 直到能达到一个较满意的精度为止。若最终对X共提取了m个成分t1, …, tm, 偏最小二乘回归将通过实施yk对t1, …, tm的回归, 然后再表达成yk关于原变量x, …, x的回归方程 (k=1, 2, …, q) 。

1.2交叉有效性检验

将所有n个样本点分成两部分:第一部分除去某个样本点i的所有样本点集合 (共含n-1个样本点) , 用这部分样本点并使用h个成分拟合一个回归方程;第二部分是把刚才被排除的样本点i代入前面拟合的回归方程, 得到yj在样本点i上的拟合值y hj (-i) 。对于每一个i=1, 2, …, n, 重复上述测试, 定义Y的预测误差平方和为PRESSh, 有PRESSh=, 显然, 如果回归方程的稳健性不好, 误差就很大, 它对样本点的变动就会十分敏感, 这种扰动误差的作用, 就会加大PRESSh的值。

1.3精度分析

在偏最小二乘回归计算过程中, 所提取的自变量成分th一方面可以尽可能多地代表X的变异信息, 另一方面又尽可能与Y相关联, 解释Y中的信息, 用r (xi, xj) 表示两个变量之间的相关关系, 定义t1, t2, …, th对各变量的解释能力。

对自变量X的累计解释能力:

对自变量Y的累计解释能力:

定义为变量投影重要性指标用来描述自变量对因变量的解释能力, 其中whj表示轴wh的第j个分量, VIPj的意义:由于xj对Y的解释是通过th来传递的, 如果th对Y的解释能力很强, 而xj在构造th时, 又起到了非常重要的作用, 则可以认为xj对Y具有相当大的解释能力。

2鱼雷管峰压预测

2.1建立模型

发射鱼雷时, 鱼雷从静止开始, 可以通过不同的运动规律, 进行舰艇鱼雷发射装置管内运动参数的设计。分别记发射气瓶的容积记为x1, 发射气瓶的初始压力为x2, 发射阀喉部的通径为x3, 文献[1]分析了内弹道试验实测数据表1。

建立了鱼雷管峰压的线性回归模型为

Pm=-0.078 3x1+0.035x2+0.048 9x3。

本文利用偏最小二乘法建立模型。首先进行交叉有效性检验, 可得PRESSh=[0.5136 0.21330.2062]', 当h=3时, PRESSh达到最小, 模型的预测能力最好, 因此取3个主成分t1, t2, t3建立模型, 得到鱼雷管峰压的偏最小二乘回归标准化变量模型:

Pm=0.153 7x1+0.795 9x2+0.367 6x3 (1)

转化为原始变量回归方程为:

Pm=-6.610 6+0.124 9x1+0.031 8x2+0.175 9x3 (2)

2.2模型评价

2.2.1累计解释能力分析

计算t1, t2, t3的累计解释能力数值见表2。可以看出, 当取3个主成分时, t1, t2, t3对x1、x2和x3的累计解释能力达到100%, 对全部自变量x的累计解释能力达到100%, 对因变量y的累计解释能力均超过98%, 均达到了很高的解释水平, 说明用偏最小二乘法拟合出的回归模型能概括原始数据所携带大部分的信息。

2.2.2变量投影重要性指标

经计算, 得变量投影重要性指标VIP=[0.1834 1.463 1 0.908 8]', 对鱼雷管峰压的影响程度依次为发射气瓶的初始压力、发射阀喉部的通径和发射气瓶的容积。

2.2.3相关关系分析

图1给出了三个主成分与u1的平面图。从图中可以看出, t1与u1存在明显的线性关系, 说明发射气瓶的容积、发射气瓶的初始压力、发射阀喉部的通径与鱼雷管峰压有显著的相关关系, 采用偏最小二乘回归方法建立的线性模型 (1) 、 (2) 是合理的, t2、t3与u1也存在着一定的线性关系, 但已经很弱了。

2.2.4 鱼雷管峰压拟合结果

鱼雷管峰压的实测值与拟合值的对比情况如表3所示。从中可以看出, 偏最小二乘回归拟合值好于文献[1]的回归值, 相对误差较小, 表明偏最小二乘回归模型能较好的反映鱼雷管峰压与气瓶压力、发射阀通径、气瓶容积的关系, 预测效果较好。

鱼雷管峰压实测值与拟合值对比图如图2所示。可以看出, 所有的样本点均排列在图中对角线的附近, 模型的拟合效果很好。

3 结论

本文采用偏最小二乘回归方法, 经交叉有效性检验, 建立了鱼雷管峰压与气瓶容积、气瓶压力和发射阀通径的关系模型。通过检验, 该模型对自变量和因变量的解释能力较高, 通过与实测值的对比, 该模型拟合精度较高, 具有较好的预测能力。

摘要：采用偏最小二乘回归方法, 经交叉有效性检验建立了鱼雷管峰压的回归方程。通过实验数据比较, 该模型拟合精度较强, 说明偏最小二乘回归是估算舰艇鱼雷发射内弹道特征参数的一种有效方法。

关键词：鱼雷管峰压,偏最小二乘,预测,误差

参考文献

[1] 李志华, 申志伟. 鱼雷内弹道特征参数估计的回归分析. 舰船电子工程, 2009 (12) :188—191

[2] 乔汝椿. 潜艇鱼雷发射装置内弹道特性初步设计. 水中兵器, 2003; (2) :51—55.

[3] 王惠文. 偏最小二乘回归方法及其应用. 北京:国防工业出版社, 1999

[4] 罗批, 郭继昌, 李锵, 等. 基于偏最小二乘回归建模的探讨. 天津大学学报, 2002;35 (6) :783—786

[5] 贾子英, 陈松辉, 黄华. 基于偏最小二乘的两栖突击车作战效能预测. 舰船电子工程, 2009; (1) :119—121

[6] 巴胜生. 基于偏最小二乘回归的军用飞机效能预测. 自动化技术与应用, 2006;25 (10) :40—42

[7] 钱筱丹, 黎放, 卞金露. 偏最小二乘回归在舰船维修费用预测中的应用. 舰船科学技术, 2007;29 (4) :98—100

[8] 任震, 张静伟, 张晋听. 基于偏最小二乘法的设备故障率计算. 电网技术, 2005;29 (5) :12—16

[9] 王文圣, 丁晶, 赵玉龙, 等. 基于偏最小一乘回归的年用电量预测研究. 中国电机工程学报, 2003;23 (10) :17—21

回归估计篇8

雇员信息表1中,每个雇员有四个属性:ID,Degree,Experience和Salary。其中有一个雇员的属性含有空值,本文将对这个空值信息进行估计。先用模糊聚类方法,将不含空值的22个雇员根据其属性进行分簇,然后将含空值的雇员划归到与他最相似簇中,最后用线性回归的方法进行计算,最终估计出空值。

因为在属性值Degree中,(Ph.D.)>(Master)>(Bachelor)。为了实验方便,我们将Ph.D.赋值为3,Master赋值为2,Bachelor赋值为1。这样的赋值对空值估计来说可以达到最佳的估计效果(Chen and Yeh 1997)。同时,将属性Salary的数值转换在0,1之间。

1 模糊聚类概述

聚类就是按照事物间的相似性进行区分和分类的过程。传统的聚类分析是一种硬划分,具有非此即彼的性质,Zadeh[7]提出的模糊理论为这种软划分提供了有力的分析工具,人们开始用模糊的方法来处理聚类问题,并称之为模糊聚类分析。模糊聚类最终将给定的样本划分成k个簇(C1,C2,…,Ck),如下:

其中1≤k≤m,Sij表示Ci中的第j个元素。

模糊聚类的实现方法主要是根据研究对象本身的属性来构造模糊矩阵,并在此基础上根据一定的隶属度来确定聚类关系,也就是用模糊数学的方法把样本之间的模糊关系定量的确定,从而客观且准确地进行聚类。这里,我们采用基于等价关系[8]的模糊聚类方法来处理,用MATLAB这一数学工具编程实现。

2 用MATLAB进行模糊聚类

MATLAB是目前国际上流行的科学计算软件,它具有强大的矩阵计算和数据可视化能力。用MATLAB模糊聚类分析过程的两大步骤:求解模糊相似矩阵和建立模糊等价矩阵。

2.1 建立模糊相似矩阵算法

用绝对值减数法[2]进行标定:,其中0≤c≤1适当选取c的值,使得rij在[0,1]中分散开来,得到模糊相似矩阵R。

基于属性值Salary,应用模糊聚类算法,将不包含空值的元组分簇。取表中的salary作为样本集合:S={S1,S2,......,S21}。

程序算法如下:

2.2 建立模糊等价矩阵

程序算法如下:

计算λ截矩阵:取λ=0.99时,算法如下:

得到相似矩阵R(略)。

2.3 聚类结果

根据矩阵得到聚类结果如表3。

此结果将表2中,把含空值以外的22个元组信息分成12个簇:(C1,C2,…,C12)。簇的中心ci可以表示成:

其中n表示Ci中元素的个数,Dij表示Ci中第j个元素Degree的属性值。Eij代表Ci中第j个元素Experience的属性值。Sij表示Ci中第j个元素Salary的属性值。

3 用线性回归求空值

用线性回归计算空值的公式如下:

是待估计的空值,Di_center,Ei_center,Si_center)是包含空值的元组被划归到的簇的中心。

βDi和βEi是属性Degree以及Experience值的单位变化对属性Salary值的偏移量,即回归系数。Wd和We是根据βDi和βEi计算的系数。

3.1 根据包含空值的元组的其它属性将其划归到最相似的簇中

用表示含空值的元组,然后计算和(Di_center,Ei_center,Si_center)之间的欧拉距离。计算公式如下:

例如含空值的与C5的中心c5的欧拉距离:Dist5=((2-1.66667)2+(0.45-0.47)2)1/2=0.333933.依此类推,算出每个Ci的Disti。如果Distp=Min(Dist1,Dist2,…,Distk),其中1≤p≤k。则将(D,E,S^)作为Ci中的元素。结果显示Dist5为Dist1,Dist2,…,Dist12中的最小值。最后将(D,E,S^)划归到簇C5中。

3.2 用回归系数来衡量其余属性值对空值属性值的影响力

计算属性Degree以及Experience值的单位变化对属性Salary值的偏移量。假设Ci中的第j个元数为(Dij,Eij,Sij),1≤j≤n,其中n代表Ci中的元素个数,则有:

如βD5和βE5计算过程和结果如下:

3.3 估计属性Salary的值线性回归公式如下:

最终估计值为0.57709×70000=40396.1

我们将提出的方法估计表中的每个元组的属性值Salary,见表4。

对空值估值算法,本文采用MAER(e^)来衡量新方法的估值准确性。定义如下:

其中为属性Salary的估计值,Sori为表中初始值,m为元组个数。

表4给出了每个估计值的误差以及和前人方法的比较结果。数据表明该方法有较低的误差率。

4 结束语

本文介绍了估计数据库空值的新方法。使用了Matlab这一数学工具编程,实现了模糊聚类分析方法。同时将得到的结果通过隶属度的加权进行了修正,得到了误差率更低的聚类结果。然后根据包含空值的元组的其它属性将其划归到最相似的簇中,最后再用线性回归法对空值进行估计。在附表4中给出了新方法与前面方法的比较结果和误差率。

参考文献

[1]苏金明,王永利MATLAB7.0实用指南[M].北京:电子工业出版社,2004:66-68.

[2]罗承忠.模糊集引论[M].北京:北京师范大学出版社,2005:105-111.

[3]Draper NR,Smith H.Applied regression analysis.Wiley,New York,NY,1998.

[4]Wang JW,Cheng CH.An efficient method for estimating null values in relational databases.Knowl Inf Syst,2006.

[5]Chen SM,Hsiao HR.A new method to esti-mate null values in relational database systems b-ased on automatic clustering techniques.Information Science.2005,169(47-69).

[6]Chen SM,Lee SW.A new method to generate fuzzy rules from relational database systems for estimating null values.Cybern Syst,2003,Int J34:33-57.

[7]Zadeh L A.Fuy sets as a basis for a theory of possibiliy.FSS,1,1978:3-28.

回归估计篇9

关键词：系统谐波阻抗,偏最小二乘回归,用户谐波发射水平

在电力系统中, 供电系统谐波的定义是对周期性非正弦电量进行傅立叶级数分解, 除了得到与电网基波频率相同的分量, 还可以得到一系列大于电网基波频率的分量, 这一部分分量即被成为谐波。谐波频率与几波频率的比值称之为谐波次数。电网中有时也存在非整数倍谐波, 称之为非谐波或者分数谐波, 谐波实际上是一种干扰量, 使电网受到污染。

1 谐波的危害和抑制

1.1 谐波的危害

理想的公用电网所提供的电压应该是单一而固定的频率以及规定的电压幅值。谐波电流和谐波电压的出现, 对公用电网是一种污染。谐波对电网的危害十分严重, 使电力系统中电能的生产、传输和利用的效率降低, 使电气设备过热、产生振动和噪声, 并使绝缘老化, 使用寿命缩短, 甚至发生故障或者烧毁等等。谐波可引起电力系统局部并联谐振或串联谐振, 使谐波含量放大, 造成电容器等设备烧毁。谐波还会引起继电保护和自动装置误动作, 使电能计量出现混乱。对于电力系统外部, 谐波对通信设备和电子设备会产生严重干扰。

1.2 谐波的抑制

对于解决电力系统中电力装备和其他谐波来源的污染问题, 有两条可行的方法, 首先是装设一种谐波补偿装置来补偿谐波, 这一方法适用于各种来源的谐波, 其次是对电力系统中的电力设备进行种种技术上的改造, 尽量使其不产生谐波, 并且始终将功率因数控制到1, 这一条方法只适用于被谐波污染的电力系统中电子设备。在传统方法中, 装设谐波补偿装置的方法就是采用LC调谐波滤波器, 这种方法技能补偿谐波, 又能补偿无功功率, 而且LC调谐波滤波器结构简单, 操作方便, 因此一直被广泛使用。而LC调谐波滤波器的主要缺点则是对谐波的补偿特性会受到电网的阻抗和电力系统的运行状态的影响, 比较容易的和电力系统发生联通谐振, 致使谐波不仅没有被补偿, 反而持续放大, 最终导致LC滤波器过载运行直至烧毁。而且LC滤波器只能补偿固定频率的谐波, 补偿效率也不高。

2 偏最小二乘回归法对谐波检测的运用

对于电力系统的运行, 谐波的治理已经是必然要面对的重点问题。谐波的定位有两种方法, 第一种是将电力系统分成供电侧和用户侧两侧, 然后根据相应的等效电路模型确定发出主谐波源存在的那一侧, 这种方法成为基于等小电路吗, 模型的定位法, 根据定位依据的不同, 可以分为功率定位法、阻抗定位法、灵敏度定位法;第二种是估计电力系统整体谐波的方法, 计算出系统中各个节点的谐波电压一直各个支路的谐波电流, 从计算结果中判断出哪条支路或者哪个节点含有谐波, 根据不同的量测量的选择, 可以分为功率量测、电压量测和电流量测三种。

谐波源的监测就是指定量的计算出系统侧与用户侧对公共连接点的谐波畸变电压的大小, 并明确估算出用户侧谐波的发射水平。国内外对于谐波发射水平的估计方法主要还是围绕着对系统和用户谐波波阻抗的估算来展开。

偏最小二乘回归法, 是一种新型的多元统计数据分析方法, 它主要研究的是多应变量对多自变量的回归剑魔, 特别当各变量内部高度线性相关时, 用偏最小二乘回归法更有效。此外, 偏最小二乘回归法还可以较好的解决了样本个数少于变两个数等问题。偏最小二乘法是集主成分分析、典型相关分析和多元线性回归分析3种分析方法的优点于一身。它与主成分分析法都试图提取出反映数据变异的最大信息, 但主成分分析法只考虑一个自变量矩阵, 而偏最小二乘法还有一个“响应”矩阵, 因此具有预测功能。研究认为, 集多元线性回归分析、典型相关分析、主因子分析等方法于一体的偏最小二乘回归方法 (PLS) 更适用于FM分析, 可以避免数据非正态分布、因子结构不确定性 (factor indeterminacy) 和模型不能识别等潜在问题。偏最最小二乘回归与传统的多元线性回归模型相比, 有以下几个突出的特点:

(1) 能够在自变量存在严重的相关性的条件下进行回归建模;

(2) 允许在样本点个数少于变量个数的条件下进行回归建模;

(3) 偏最小二乘回归在最终模型中将包含原有的所有自变量;

(4) 偏最小二乘回归模型更易于辩识系统信息与噪声 (甚至一些非随机的噪声) ;

(5) 在偏最小二乘回归模型中, 每一个自变量的回归系数将更容易解释。

目前, 偏最小二乘回归计算方法在国内已经开始逐步实施, 该方法应用于负荷模型预测的研究成果, 并且在过去“二元回归”的方法的基础上, 已有将偏最最小二乘回归方法初步应用于配电网谐波源定位与检测中, 充分的发挥算法本身的特点, 对电力系统数据建模有了新的思路。在应用偏最小二乘回归方法进行谐波源的定位时, 首先应用PMU (相量测量单元) 对电力系统等效电路各个节点的电压和电流相量进行同步测量, 使得测量结果更加准确, 从而给偏最小二乘回归法的应用奠定了更好的基础, 使得计算结果更加符合实际情况, 以使谐波源的定位更加准确。谐波阻抗与谐波发射水平评估是近年来国内外较为关注的谐波源检测方面的研究内容。利用偏最小二乘回归方法可得到较好的评估效果。采用统计回归和其它建模方法, 可在大量相关的实验数据基础上对复杂系统进行建模, 但普通的多元统计回归方法在实际系统建模中往往存在一些难以解决的问题。传统的最小二乘法在其参数估计式的计算中, 要求其计算公式必须是可逆矩阵, 估计式才有意义。所以当模型变量的多重相关性严重时, 或者当系统中样本容量少于变量个数时, 参数估计一般就会失效。采用人工神经网络建模时, 网络模型的选取通常只能依据经验或采用随机试探的方法, 具有一定的随意性;另外, 由于网络采用“黑箱”结构, 它对所描述对象的输入输出变量之间的关系往往缺乏很好的解释性。这些问题在一定程度上制约了神经网络在系统建模方面的进一步发展。偏最小二乘回归方法借助提取主元的思路, 有效地提取对系统解释性最强的综合信息, 从而实现对高维数据空间的降维处理, 较好地克服变量多重相关性在系统建模过程中的不良影响。二元回归法是一种估算电力系统谐波阻抗和用户谐波发射水平的方法, 但在二元回归方程中, 两个自变量的相关系数较大, 存在严重的多重相关性, 这样会使最小二乘法失效, 得到的回归模型的拟合效果不好。为了较好地解决多重相关性问题, 基于电网各参数的复数关系的二元回归方程, 可以用偏最小二乘法来求解系统阻抗和用户谐波发射水平。综上所述, 根据偏最小二乘的基本原理、建模基本思想, 以及该方法在配电网谐波电源定位与检测中的应用, 已经可以证实偏最小二乘法在配电网谐波源定位与检测中得到了较好的效果。

由于实习系统中谐波源发射水平较低, 在数据处理的精度上出现误差相对较大的情况, 当系统中谐波含量越大, 偏最小二乘回归模型估算方法在数据建模上的效果将更加明显。

3 结语

随着电力企业的不断发展壮大, 为了适应社会和人民生活工作的要求, 必然要对电能的质量做出更高的要求, 谐波治理是电力系统中影响电能质量的重点, 应用偏最小二乘回归模型可以较好的估算出系统谐波阻抗以及用户谐波发射水平, 对于分清楚电力系统中系统侧和用户侧的对于电能质量恶化的责任。与传统方法相比, 偏最小二乘回归法算法集多元性回归、典型相关分析和主成分分析的基本功能为一体, 可以同时实现回归建模、数据结构简化以及变量间相关性分析, 能够有效的解决谐波检测中变量相关性给系统建模带来的误差, 并且更加的统一辨认出系统信息与造成, 对于电力系统中谐波的检测和定位有很大的帮助。

参考文献

[1]黄舜, 徐永海.基于偏最小而成回归的系统谐波阻抗与谐波发射水平的评估方法[J].中国电机工程学报, 2007, 27 (1) :94-97.

[2]侯丽丽.配电网谐波源的定位[D].山东大学, 2009.

【回归估计】推荐阅读：

线性回归问题中的一个估计08-09

密度估计05-29

背景估计06-03

稀疏估计06-15

精度估计07-12

转速估计07-29

优化估计08-04

安全估计08-05