天文算法2

May 13, 2018 5 分钟阅读时间 golang, 算法

第四章拟合 Curve Fitting

1. 什么是拟合

曲线拟合(Curve Fitting)的数学定义是指用连续曲线近似地刻画或比拟平面上一组离散点所表示的坐标之间的函数关系，是一种用解析表达式逼近离散数据的方法。曲线拟合通俗的说法就是“拉曲线”，也就是将现有数据透过数学方法来代入一条数学方程式的表示方法。科学和工程遇到的很多问题，往往只能通过诸如采样、实验等方法获得若干离散的数据，根据这些数据，如果能够找到一个连续的函数（也就是曲线）或者更加密集的离散方程，使得实验数据与方程的曲线能够在最大程度上近似吻合，就可以根据曲线方程对数据进行数学计算，对实验结果进行理论分析，甚至对某些不具备测量条件的位置的结果进行估算。

2. 最小二乘拟合

　　考虑 N 个数据点，它们的坐标是$(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)$ 假设这些值中的 X 是严格的精确值，Y 的值是测量值(含有一些误差)。对于一个给定的 X，如$x_1$，对应的值$y_1$与曲线 C 上对应的 Y 值将存在一个差值$d_1$，我们称这个差值为偏差、误差或残差，它可能是正、负或零。类似的，$x_2,\dots,x_n$,对应的差值为 $d_2,\dots,d_n$。

　　我们用 $d_1^2+d_2^2+\cdots+d_n^2$ 作为衡量曲线 C 拟合的“最佳”程度，这个值越小越好，越大则越不好。因此，我们做以下定义：任何一种类型的曲线，它们都有一个共同的特性，当 $$\sum_{i=1}^n d_i^2$$最小时，称为最佳拟合曲线。

　　一条曲线具有这一特性时，称之为“最小二乘拟合”，这样的曲线称为“最小二乘曲线”。

2.1 线性拟合

线性方程

假设拟合的曲线为直线$y=ax+b$,则最小二乘差方和为： $$\sum_{i=1}^n[y_i-(ax_i+b)]^2$$ 所以问题转化为求解上述关于a,b的二元函数的最小值

对上式中a求偏导

\begin{multline} \shoveleft \begin{aligned} \frac \partial {\partial a}\sum_{i=1}^n[y_i-(ax_i+b)]^2 &= -2\sum_{i=1}^n({x_iy_i}-a{x_i^2}-b {x_i})\\[2ex] & =-2\sum_{i=1}^n{x_i}[{y_i}-(a{x_i}+b)]\\[2ex] \end{aligned} \end{multline}

\begin{multline} \shoveleft \text {由极值的必要条件可知，}\\[2ex] \shoveleft -2\sum_{i=1}^n{x_i}[{y_i}-(a{x_i}+b)] = 0\\[2ex] \shoveleft \therefore \sum_{i=1}^n{x_iy_i}-a\sum_{i=1}^n{x_i^2}-b\sum_{i=1}^n{x_i}=0 \end{multline}

同理对b求偏导可得： \begin{align} \sum_{i=1}^n{y_i}-a\sum_{i=1}^n{x_i}-bn=0 \end{align}

求解关于a,b 的二元一次方程组， \begin{cases} \sum\limits_{i=1}^n{x_iy_i}-a\sum\limits_{i=1}^n{x_i^2}-b\sum\limits_{i=1}^n{x_i}=0 \\[2ex] \sum\limits_{i=1}^n{y_i}-a\sum\limits_{i=1}^n{x_i}-bn=0 \end{cases}

得 \begin{align} a &= \frac {\sum\limits_{i=1}^n{x_i}\sum\limits_{i=1}^n{y_i}-n\sum\limits_{i=1}^n{x_iy_i}}{ (\sum\limits_{i=1}^n{x_i})^2-n\sum\limits_{i=1}^n{x_i}^2}\\[2ex] b &= \frac {\sum\limits_{i=1}^n{x_i}\sum\limits_{i=1}^n{x_iy_i}-\sum\limits_{i=1}^n{y_i} \sum\limits_{i=1}^n{x_i}^2}{(\sum\limits_{i=1}^n{x_i})^2-n\sum\limits_{i=1}^n{x_i}^2} \end{align}
皮尔逊相关系数r

定义

由定义可知， \begin{align} r=\frac {\sum\limits_{i=1}^n(x_i-\bar x)(y_i-\bar y)}{\sqrt {\sum\limits_{i=1}^n(x_i-\bar x)^2}\sqrt {\sum\limits_{i=1}^n(y_i-\bar y)^2}}, \\[2ex] \text {其中 } \bar x=\frac 1n\sum_{i=1}^nx_i,\bar y=\frac 1n\sum_{i=1}^ny_i \end{align}

\begin{multline} \shoveleft \begin{aligned} 化简 \ \ & \sum\limits_{i=1}^n(x_i-\bar x)(y_i-\bar y)\\[2ex] = & \sum\limits_{i=1}^n(x_iy_i-x_i\bar y-\bar xy_i+\bar x\bar y)\\[2ex] = & \sum\limits_{i=1}^nx_iy_i - \sum\limits_{i=1}^nx_i\bar y - \sum\limits_{i=1}^n\bar xy_i+\sum\limits_{i=1}^n\bar x\bar y\\[2ex] = & \sum\limits_{i=1}^nx_iy_i - \frac 1n\sum\limits_{i=1}^ny_i\sum\limits_{i=1}^nx_i - \frac 1n\sum\limits_{i=1}^nx_i\sum\limits_{i=1}^ny_i + \frac 1n\sum\limits_{i=1}^nx_i\sum\limits_{i=1}^ny_i\\[2ex] = & \sum\limits_{i=1}^nx_iy_i - \frac 1n\sum\limits_{i=1}^ny_i\sum\limits_{i=1}^nx_i\\[2ex] 化简 \ \ & \sum\limits_{i=1}^n(x_i-\bar x)^2\\[2ex] = & \sum\limits_{i=1}^n(x_i^2-2x_i\bar x+\bar x^2)\\[2ex] = & \sum\limits_{i=1}^nx_i^2 - 2\bar x\sum\limits_{i=1}^nx_i + n\bar x^2\\[2ex] = & \sum\limits_{i=1}^nx_i^2 - 2n\bar x^2 + n\bar x^2\\[2ex] = & \sum\limits_{i=1}^nx_i^2 - n\bar x^2\\[2ex] = & \sum\limits_{i=1}^nx_i^2 - \frac 1n(\sum\limits_{i=1}^nx_i)^2\\[2ex] 同上 \ \ & \sum\limits_{i=1}^n(y_i-\bar y)^2\\[2ex] = & \sum\limits_{i=1}^ny_i^2 - \frac 1n(\sum\limits_{i=1}^ny_i)^2\\[2ex] \end{aligned} \end{multline} \begin{multline} \shoveleft \begin{aligned} \therefore \ \ r&=\frac {\sum\limits_{i=1}^nx_iy_i - \frac 1n\sum\limits_{i=1}^ny_i\sum\limits_{i=1}^nx_i}{\sqrt{\sum\limits_{i=1}^nx_i^2 - \frac 1n(\sum\limits_{i=1}^nx_i)^2}\sqrt{\sum\limits_{i=1}^ny_i^2 - \frac 1n(\sum\limits_{i=1}^ny_i)^2}}\\[2ex] & =\frac {n\sum\limits_{i=1}^nx_iy_i - \sum\limits_{i=1}^ny_i\sum\limits_{i=1}^nx_i}{\sqrt{n\sum\limits_{i=1}^nx_i^2 - (\sum\limits_{i=1}^nx_i)^2}\sqrt{n\sum\limits_{i=1}^ny_i^2 - (\sum\limits_{i=1}^ny_i)^2}} \end{aligned} \end{multline}

　　这个系数介于+1 到-1 之间。如果值为+1 或-1，说明 x 和 y 之间有完全的线性关系，所有的点(x,y)精确的在同一条直线上。如果 r = +1，y 随 x 单调递增，如果 r = -1，y 随 x 单调递减。

代码

　　到这里我们就已经推导出求解线性方程和相关系数的公式，再结合代码看看

// Linear fits a line to sample data.
//
// Argument p is a list of data points.  Results a and b are coefficients
// of the best fit line y = ax + b.
// 求解线性拟合直线
// sx = ∑x sy = ∑y sxy = ∑xy sx2 = ∑x^2
func Linear(p []struct{ X, Y float64 }) (a, b float64) {
    var sx, sy, sx2, sxy float64
    for i := range p {
        x := p[i].X
        y := p[i].Y
        sx += x
        sy += y
        sx2 += x * x
        sxy += x * y
    }
    n := float64(len(p))
    d := n*sx2 - sx*sx
    // (4.2) p. 36
    a = (n*sxy - sx*sy) / d
    b = (sy*sx2 - sx*sxy) / d
    return
}

// CorrelationCoefficient returns a correlation coefficient for sample data.
// 求解相关系数 r
func CorrelationCoefficient(p []struct{ X, Y float64 }) float64 {
    var sx, sy, sx2, sy2, sxy float64
    for i := range p {
        x := p[i].X
        y := p[i].Y
        sx += x
        sy += y
        sx2 += x * x
        sy2 += y * y
        sxy += x * y
    }
    n := float64(len(p))
    // (4.3) p. 38
    return (n*sxy - sx*sy) / (math.Sqrt(n*sx2-sx*sx) * math.Sqrt(n*sy2-sy*sy))
}

2.2 二次曲线拟合

二次方程$y = ax^2 + bx + c$

　　假设我们希望画一条逼近 N 个点的最佳二次曲线：$y = ax^2 + bx + c$ 这是一个纵轴的抛物线。同一次直线类似，差方和为$$\sum_{i=1}^n[y_i-(ax_i^2+bx_i+c)]^2$$ 依次对上式中的a,b,c求偏导，得 \begin{align} \sum_{i=1}^n(-2x_i^2y_i+2ax_i^4+2bx_i^3+2cx_i^2) = 0\\[2ex] \sum_{i=1}^n(-2x_iy_i+2ax_i^3+2bx_i^2+2cx_i) = 0\\[2ex] \sum_{i=1}^n(-2y_i+2ax_i^2+2bx_i+2c) = 0\\[2ex] \end{align} \begin{multline} \begin{aligned} \therefore & \sum_{i=1}^nx_i^2y_i=\sum_{i=1}^nax_i^4+\sum_{i=1}^nbx_i^3+\sum_{i=1}^ncx_i^2\\[2ex] & \sum_{i=1}^nx_iy_i=\sum_{i=1}^nax_i^3+\sum_{i=1}^nbx_i^2+\sum_{i=1}^ncx_i\\[2ex] & \sum_{i=1}^ny_i=\sum_{i=1}^nax_i^2+\sum_{i=1}^nbx_i+nc\\[2ex] \end{aligned} \end{multline}

求解以下三元一次方程组 \begin{cases} \sum\limits_{i=1}^nx_i^2y_i=\sum\limits_{i=1}^nax_i^4+\sum\limits_{i=1}^nbx_i^3+\sum\limits_{i=1}^ncx_i^2\\[2ex] \sum\limits_{i=1}^nx_iy_i=\sum\limits_{i=1}^nax_i^3+\sum\limits_{i=1}^nbx_i^2+\sum\limits_{i=1}^ncx_i\\[2ex] \sum\limits_{i=1}^ny_i=\sum\limits_{i=1}^nax_i^2+\sum\limits_{i=1}^nbx_i+nc\\[2ex] \end{cases}

得 \begin{multline} \shoveleft \begin{aligned} a &= \frac {NQV+PRT+PQU-Q^2T-P^2V-NRU}{D}\\[2ex] b &= \frac {NSU+PQV+QRT-Q^2U-PST-NRV}{D}\\[2ex] c &= \frac {QST+QRU+PRV-Q^2V-PSU-R^2T}{D}\\[2ex] \end{aligned} \end{multline}

\begin{multline} \shoveleft \begin{aligned} 其中\ &P =\sum\limits_{i=1}^nx_i,\ Q=\sum\limits_{i=1}^nx_i^2,\ R=\sum\limits_{i=1}^nx_i^3,\ S=\sum\limits_{i=1}^nx_i^4,\\[2ex] &T=\sum\limits_{i=1}^ny_i,\ U=\sum\limits_{i=1}^nx_iy_i,\ V=\sum\limits_{i=1}^nx_i^2y_i,\\[2ex] &D = NQS+2PQR-Q^3-P^2S-NR^2 \end{aligned} \end{multline}

　　至此我们的程序就能很方便的通过上述公式，对这 N 个点进行一次完整的遍历，计算上述 P,Q,R,S,T,U,V,D,即可求得拟合曲线方程。

代码

// Quadratic fits y = ax² + bx + c to sample data.
//
// Argument p is a list of data points.  Results a, b, and c are coefficients
// of the best fit quadratic y = ax² + bx + c.
// 求解二次拟合曲线系数
func Quadratic(p []struct{ X, Y float64 }) (a, b, c float64) {
    var P, Q, R, S, T, U, V float64
    for i := range p {
        x := p[i].X
        y := p[i].Y
        x2 := x * x
        P += x
        Q += x2
        R += x * x2
        S += x2 * x2
        T += y
        U += x * y
        V += x2 * y
    }
    N := float64(len(p))
    // (4.5) p. 43
    D := N*Q*S + 2*P*Q*R - Q*Q*Q - P*P*S - N*R*R
    // (4.6) p. 43
    a = (N*Q*V + P*R*T + P*Q*U - Q*Q*T - P*P*V - N*R*U) / D
    b = (N*S*U + P*Q*V + Q*R*T - Q*Q*U - P*S*T - N*R*V) / D
    c = (Q*S*T + Q*R*U + P*R*V - Q*Q*V - P*S*U - R*R*T) / D
    return
}

一般曲线拟合（多重线回归）

　　最佳线性拟合的原理可以被扩展到其它函数，这个函数可以含有超过两个未知的线性系数。让我们考虑三个函数的线性组合的情况。假设我们已知：$$y = af_0(x) + bf_1(x) + cf_2(x)$$ 式中$f_0$、$f_1$ 和 $f_2$ 是三个关于 x 的已知函数，但系数 a、 b 和 c 是未知的。此外，假设已知 3 个 x 对应的 y 值。那么系数 a、b、c 可按如下得到。求和计算： \begin{align} M &= \sum f_0^2 &U &= \sum yf_0\\[2ex] P &= \sum f_0f_1 &V &= \sum yf_1\\[2ex] Q &= \sum f_0f_2 &W &= \sum yf_2\\[2ex] R &= \sum f_1^2 \\[2ex] S &= \sum f_1f_2 \\[2ex] T &= \sum f_2^2 \\[2ex] \end{align} $$D = MRT+2PQS-MS^2-RQ^2-TP^2$$ 那么： \begin{align} a &= \frac {U(RT-S^2)+V(QS-PT)+W(PS-QR)}{D}\\[2ex] b &= \frac {U(SQ-PT)+V(MT-Q^2)+W(PQ-MS)}{D}\\[2ex] c &= \frac {U(PS-RQ)+V(PQ-MS)+W(MR-P^2)}{D} \end{align} 另一种特殊情况，考虑 y=af(x)，只有一个未知系数。我们容易得到： $$a = \frac {\sum yf}{\sum f^2}$$

代码

// Func3 implements multiple linear regression for a linear combination
// of three functions.
//
// Given sample data and three functions in x, Func3 returns coefficients
// a, b, and c fitting y = aƒ₀(x) + bƒ₁(x) + cƒ₂(x) to sample data.
// 多重线性回归
func Func3(p []struct{ X, Y float64 }, f0, f1, f2 func(float64) float64) (a, b, c float64)  {
    var M, P, Q, R, S, T, U, V, W float64
    for i := range p {
        x := p[i].X
        y := p[i].Y
        y0 := f0(x)
        y1 := f1(x)
        y2 := f2(x)
        M += y0 * y0
        P += y0 * y1
        Q += y0 * y2
        R += y1 * y1
        S += y1 * y2
        T += y2 * y2
        U += y * y0
        V += y * y1
        W += y * y2
    }
    // (4.7) p. 44
    D := M*R*T + 2*P*Q*S - M*S*S - R*Q*Q - T*P*P
    a = (U*(R*T-S*S) + V*(Q*S-P*T) + W*(P*S-Q*R)) / D
    b = (U*(S*Q-P*T) + V*(M*T-Q*Q) + W*(P*Q-M*S)) / D
    c = (U*(P*S-R*Q) + V*(P*Q-M*S) + W*(M*R-P*P)) / D
    return
}

// Func1 fits a linear multiple of a function to sample data.
//
// Given sample data and a function in x, Func1 returns coefficient
// a fitting y = aƒ(x).
func Func1(p []struct{ X, Y float64 }, f func(float64) float64) float64 {
    var syf, sf2 float64
    // (4.8) p. 45
    for i := range p {
        f := f(p[i].X)
        y := p[i].Y
        syf += y * f
        sf2 += f * f
    }
    return syf / sf2
}

天文算法

最近更新于 May 13, 2018