@betasy 2015-03-19T04:47:05.000000Z 字数 4203 阅读 6206

多元线性回归求解任意多边形的最优椭圆（Python实现）

算法

1. INPUT 多边形边界序列数组

现在使用的数据，是夜间灯光数据边界查询后得到的某区域对象多边形的边界像元坐标点对。形如：
B=[ [1,2,3,4,5,6], [3,2,4,5,4,3] ]
的二维数组。

2. 加载数组到内存

使用numpy模块提供的loadtxt的方法加载数组，代码从略。

3. 建立任意椭圆代数解析式

一般地，任意椭圆的方程可以表示为：

A x 2 + B x y + C y 2 + D x + E y + F = 0

$Ax^2+Bxy+Cy^2+Dx+Ey+F=0$
其中

A,B,C,D,E,F $A,B,C ,D, E, F$ 为参数，对于任意点

P(xi,yi) $P(x_i,y_i)$ ，相对于椭圆的位置都有三种情况：
1.在椭圆上，

Ax2i+Bxiyi+Cy2i+Dxi+Eyi+F=0 $Ax_i^2+Bx_iy_i+Cy_i^2+Dx_i+Ey_i+F=0$
2.在椭圆内，

Ax2i+Bxiyi+Cy2i+Dxi+Eyi+F<0 $Ax_i^2+Bx_iy_i+Cy_i^2+Dx_i+Ey_i+F<0$
3.在椭圆外，

Ax2i+Bxiyi+Cy2i+Dxi+Eyi+F>0 $Ax_i^2+Bx_iy_i+Cy_i^2+Dx_i+Ey_i+F>0$
设

ρ2=Ax2+Bxy+Cy2+Dx+Ey+F(1) $ρ^2=Ax^2+Bxy+Cy^2+Dx+Ey+F\qquad (1)$

4. 建立边界点与椭圆的关系

本文所求解的最优拟合椭圆需要满足以下几个条件：
1.满足椭圆定义，即 $B^2-4AC≤0$
2. $\sum_{i=1}^n (Ax_i^2+Bx_iy_i+Cy_i^2+Dx_i+Ey_i+F-ρ^2)^2$ 达到最小

5. 多元线性回归求参数

令 $Q=\sum_{i=1}^n (Ax_i^2+Bx_iy_i+Cy_i^2+Dx_i+Ey_i+F-ρ^2)^2$
要使用多元线性回归的方法拟合椭圆的参数，就要对边界点与椭圆关系的方程进行变化，以适应要求。
首先回顾一下多元线性回归模型的建立过程：
假设某一因变量 $y$ 受 $k$ 个自变量 $x_1,x_2,...,x_k$ 的影响，其 $n$ 组观测值为 $(y_a,x_{1a},x_{2a},...x_{ka}),a=1,2,...,n$ 。那么，多元线性回归模型的结构形式为：

y a = β 0 + β 1 x 1 a + β 2 x 2 a + . . . + β k x k a + ε a

$y_a=β_0+β_1x_{1a}+β_2x_{2a}+...+β_kx_{ka}+ε_a$
式中：

β0,β1,...,βk $β_0,β_1,...,β_k$ 为待定参数；

εa $ε_a$ 为随机变量。
如果

b0,b1,...,bk $b_0,b_1,...,b_k$ 分别为

β0,β1,...,βk $β_0,β_1,...,β_k$ 的拟合值，则回归方程为

y\^=b0+b1x1+b2x2+...+bkxk

$y^\^ =b_0+b_1x_1+b_2x_2+...+b_kx_k$
式中：

b0 $b_0$ 为常数；

b0,b1,b2,...,bk $b_0,b_1,b_2,...,b_k$ 称为偏回归系数。
偏回归系数

Bi(i=1,2,...,k) $B_i(i=1,2,...,k)$ 的意义是，当其它自变量

xj(j≠i) $x_j(j\not=i)$ 都固定时，自变量

xi $x_i$ 每变化一个单位而使因变量

y $y$ 平均改变的数值。
根据最小二乘法原理，

βii=0,1,2,...,k) $β_ii=0,1,2,...,k)$ 的估计值

bi(i=0,1,2,...,k) $b_i(i=0,1,2,...,k)$ 应该使

Q=∑a=1n(ya−y\^a)2=∑a=1n[ya−(b0+b1x1a+b2x2a+,...,bkxka)]2→min

$Q=\sum\limits_{a=1}^n(y_a-y^\^_a)^2=\sum\limits_{a=1}^n[y_a-(b_0+b_1x_{1a}+b_2x_{2a}+,...,b_kx_{ka})]^2\rightarrow min$
由求极值的必要条件得

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪∂Q∂b0=−2∑a=1n(ya−y\^a)=0∂Q∂bj=−2∑a=1n(ya−y\^a)xja=0(j=1,2,...,k)

$\begin{cases} \frac{\partial Q}{\partial b_0} =-2\sum\limits_{a=1}^n(y_a-y^\^_a)=0 \\ \frac{\partial Q}{\partial b_j}=-2\sum\limits_{a=1}^n(y_a-y^\^_a)x_{ja}=0\qquad (j=1,2,...,k)\ \end{cases}$
将方程组展开整理可得正规方程组，这里不再写出具体形式。
所以我们可以将

(1) $(1)$ 式看做一个观测值的通式，

∀p=(xi,yi) $\forall p = (x_i,y_i)$ 边界点,都可以建立形如

(1) $(1)$ 式的方程：

ρ 2 i = F + A x 2 i + B x i y i + C y 2 i + D x i + E y i (i = 1, 2, . . ., n)

$ρ_i^2=F+Ax_i^2+Bx_iy_i+Cy_i^2+Dx_i+Ey_i\qquad (i=1,2,...,n)$
这里可以看做

ρ2 $ρ^2$ 受5个自变量影响，则对应于多元线性回归模型的结构形式，

ρ 2 x 2 i x i y i y 2 i x i y i a = = = = = = = y a x 1 a x 2 a x 3 a x 4 a x 5 a (1, 2, 3, . . ., n)

$\begin{array}{lcl} ρ^2 & = & y_a \\ x_i^2 & = & x_{1a} \\ x_iy_i& = & x_{2a} \\ y_i^2& = & x_{3a} \\ x_i &=& x_{4a} \\ y_i &=&x_{5a}\\ a&=&(1,2,3,...,n) \end{array}$
所以引入以下向量和矩阵：

b = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ F A B C D E ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, Y = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ρ 21 ρ 22 . . . ρ 2 n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ X = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 111 ⋮ 1 x 11 x 12 x 13 ⋮ x 1 n x 21 x 22 x 23 ⋮ x 2 n x 31 x 32 x 33 ⋮ x 3 n x 41 x 42 x 43 ⋮ x 4 n x 51 x 52 x 53 ⋮ x 5 n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\boldsymbol{b}= \begin{pmatrix} F \\ A \\ B\\ C\\ D\\ E\\ \end{pmatrix}, \boldsymbol{Y}= \begin{pmatrix} ρ_1^2\\ ρ_2^2\\ .\\ .\\ .\\ ρ_n^2 \end{pmatrix} \boldsymbol{X}= \begin{pmatrix} 1 & x_{11} & x_{21}& x_{31} &x_{41}&x_{51} \\ 1 & x_{12} & x_{22}& x_{32} &x_{42}&x_{52}\\ 1 & x_{13} & x_{23}& x_{33} &x_{43}&x_{53}\\ \vdots & \vdots&\vdots&\vdots&\vdots&\vdots\\ 1 & x_{1n} & x_{2n}& x_{3n} &x_{4n}&x_{5n} \end{pmatrix}$

A = X T X = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 x 11 x 21 x 31 x 41 x 51 1 x 12 x 22 x 32 x 42 x 52 1 x 13 x 23 x 33 x 43 x 53 \dots \dots \dots \dots \dots \dots 1 x 1 n x 2 n x 3 n x 4 n x 5 n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\boldsymbol {A}=\boldsymbol{X^TX}= \begin{pmatrix} 1&1&1&\cdots&1\\ x_{11}&x_{12}&x_{13}&\cdots&x_{1n}\\ x_{21}&x_{22}&x_{23}&\cdots&x_{2n}\\ x_{31}&x_{32}&x_{33}&\cdots&x_{3n}\\ x_{41}&x_{42}&x_{43}&\cdots&x_{4n}\\ x_{51}&x_{52}&x_{53}&\cdots&x_{5n}\\ \end{pmatrix}$

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 111 ⋮ 1 x 11 x 12 x 13 ⋮ x 1 n x 21 x 22 x 23 ⋮ x 2 n x 31 x 32 x 33 ⋮ x 3 n x 41 x 42 x 43 ⋮ x 4 n x 51 x 52 x 53 ⋮ x 5 n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\begin{pmatrix} 1 & x_{11} & x_{21}& x_{31} &x_{41}&x_{51} \\ 1 & x_{12} & x_{22}& x_{32} &x_{42}&x_{52}\\ 1 & x_{13} & x_{23}& x_{33} &x_{43}&x_{53}\\ \vdots & \vdots&\vdots&\vdots&\vdots&\vdots\\ 1 & x_{1n} & x_{2n}& x_{3n} &x_{4n}&x_{5n} \end{pmatrix}$

B = X T Y = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 x 11 x 21 x 31 x 41 x 51 1 x 12 x 22 x 32 x 42 x 52 1 x 13 x 23 x 33 x 43 x 53 \dots \dots \dots \dots \dots \dots 1 x 1 n x 2 n x 3 n x 4 n x 5 n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ρ 21 ρ 22 . . . ρ 2 n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\boldsymbol{B}=\boldsymbol{X^TY}= \begin{pmatrix} 1&1&1&\cdots&1\\ x_{11}&x_{12}&x_{13}&\cdots&x_{1n}\\ x_{21}&x_{22}&x_{23}&\cdots&x_{2n}\\ x_{31}&x_{32}&x_{33}&\cdots&x_{3n}\\ x_{41}&x_{42}&x_{43}&\cdots&x_{4n}\\ x_{51}&x_{52}&x_{53}&\cdots&x_{5n}\\ \end{pmatrix} \begin{pmatrix} ρ_1^2\\ ρ_2^2\\ .\\ .\\ .\\ ρ_n^2 \end{pmatrix}$
前面说到正规方程组，正规方程组可以进一步写成矩阵形式：

A b = B (2)

$\boldsymbol{Ab=B}\qquad(2)$
求解

(2) $(2)$ 可得:

b = A - 1 B = (X T X) - 1 X T Y (3)

$\boldsymbol{b=A^{-1}}B=(X^TX)^{-1}X^TY\qquad(3)$
因此解出

(3) $(3)$ 式即可求得椭圆参数。