回归分析

回归简单来说就是由结果得到原因，是一种归纳的思想，当看到大量的事实所呈现的样态，推断出原因是如何的；当看到大量数字对是某种样态时，推断出它们之间蕴含的关系是如何的。

1. 最小二乘法

线性回归利用回归分析来确定两种或两种以上变量之间相互依赖的定量关系的一种统计分析方法。表达式：

$$y=ax+b+e$$
$e$为服从均值为0的正太分布。注意：它不是一个定值，它和$y,x$对应着出现。

示例：

由上面的打点计时计算重力加速度的实验可以得到一组数据：

其中$v_n$表示第$n$个$0.02s$的瞬时速度，是可以通过测量位移差求出来。

代码：

#coding=utf-8
import numpy as np
import matplotlib.pyplot as plt

def fun():
    # 原始数据
    x = [1, 2, 3, 4, 5, 6, 7, 8, 9]
    y = [0.199, 0.389, 0.580, 0.783, 0.980, 1.177, 1.380, 1.575, 1.771]

    A = np.vstack([x, np.ones(len(x))]).T # 垂直方向合并数组
    #A:
    #[[ 1.  1.]
    # [ 2.  1.]
    # [ 3.  1.]
    # [ 4.  1.]
    # [ 5.  1.]
    # [ 6.  1.]
    # [ 7.  1.]
    # [ 8.  1.]
    # [ 9.  1.]]

    #调用最小二乘法函数
    a, b = np.linalg.lstsq(A, y)[0]

    #转换成numpy arrary
    x = np.array(x)
    y = np.array(y)

    #画图
    plt.plot(x, y, 'o', label='Original data', markersize=10)
    plt.plot(x, a * x + b, 'r', label='Fitted line')
    plt.show()
    plt.savefig('test.png')

if __name__ == '__main__':
    fun()

运行结果：

注意

np.vstack([x, np.ones(len(x))]).T的用法：使用numpy下的vstack（垂直方向）和hstack（水平方向）函数合并数组。

>>> a = np.ones((2,2))  
>>> b = np.eye(2)  
>>> print np.vstack((a,b))  
[[ 1.  1.]  
 [ 1.  1.]  
 [ 1.  0.]  
 [ 0.  1.]]  
>>> print np.hstack((a,b))  
[[ 1.  1.  1.  0.]  
 [ 1.  1.  0.  1.]]  

这里而且是深拷贝：

>>> c = np.hstack((a,b))  
>>> print c  
[[ 1.  1.  1.  0.]  
 [ 1.  1.  0.  1.]]  
>>> a[1,1] = 5  
>>> b[1,1] = 5  
>>> print c  
[[ 1.  1.  1.  0.]  
 [ 1.  1.  0.  1.]]  

2. 残差分析

对平面上的一系列点用一条光滑的曲线连接起来的过程叫拟合，而多种拟合方法究竟哪一种画法最科学？这就要用到残差分析。也就是说在拟合的过程中，每一个点所产生的误差$e$大部分应该在0附近，而远离0误差的$e$越少越好。这就是为什么说$e$是要服从均值为0的正态分布。

假设有多个$x$和$y$的样本值，同时尝试用$y=ax+b+e$来拟合，可以得到：$|e|=|ax+b-y|$，试着把所有的$|e|$都求和。构造一个函数;

$$Q=\sum_{i=1}^n(ax_i+b-y_i)^2$$
现在就把问题转换为优化问题了，用 $Q$分别对$a$和 $b$求偏微分，满足每个偏微分方程为0的$a, b$变量的值就是要找的值。
这里用的就是残差分析。

#coding=utf-8
import numpy as np
import matplotlib.pyplot as plt


def fun():
    # 原始数据
    x = [1, 2, 3, 4, 5, 6, 7, 8, 9]
    y = [0.199, 0.389, 0.580, 0.783, 0.980, 1.177, 1.380, 1.575, 1.771]

    t1 = t2 = t3 = t4 = 0
    n = len(x)
    for i in range(n):
        t1 += y[i]
        t2 += x[i]
        t3 += x[i] * y[i]
        t4 += x[i] ** 2

    #化简后的最终形式的a,b
    a = (t1 * t2 / n - t3) / (t2 * t2 / n - t4) 
    b = (t1 - a * t2) / n

    #转换成numpy arrary
    x = np.array(x)
    y = np.array(y)

    #画图
    plt.plot(x, y, 'o', label='Original data', markersize=10)
    plt.plot(x, a * x + b, 'r', label='Fitted line')
    plt.show()
    plt.savefig('test.png')


if __name__ == '__main__':
    fun()

运行结果：

3. 非线性最小二乘法

非线性回归的情况太复杂，在上产实践中也要尽量避免使用这种模型，好在分类算法有很多，而且更多的是为了处理半结构化数据。
非线性回归一般分为一元非线性回归和多元非线性回归。

通常的做法是把非线性的回归模型经过适当的数学变换转化为线性的多元回归模型来处理。但对于另外一些非线性模型，仅仅做变量变换根本无济于事，属于前一种情况的非线性回归模型一般称为内蕴的线性回归，而后者称为内蕴的非线性回归。

示例：

根据马尔萨斯模型：

$$s=a \cdot e^{\beta t}$$
其中 $s$是人口，$t$是年份， $e$是自然常数，试推到到2030年世界人口的数量，对此分线性的多元模型可以通过变换得到如下的式子：
$$lns=\beta t+lna$$
这里就是线性模型了

代码：

#coding=utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


def fun():
    T = [1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968]
    S = [29.72, 30.61, 31.51, 32.13, 32.34, 32.85, 33.56, 34.20, 34.83]

    xdata = np.array(T)
    ydata = np.log(np.array(S)) # 这里是以e为底
    print 'xdata:',xdata
    print 'ydata:',ydata

    def func(x, a, b):
        return a + b * x

    # 使用非线性最小二乘法拟合函数
    popt, pcov = curve_fit(func, xdata, ydata) # 这里的popt是系数
    print 'popt, pcov:',popt, '\n....',pcov

    # 绘图
    plt.figure()
    plt.plot(xdata, ydata, 'ko', label="Original Noised Data")
    plt.plot(xdata, func(xdata, *popt), 'r', label="Fitted Curve")
    plt.show()
    plt.savefig('test.png')


if __name__ == '__main__':
    fun()

运行结果：

这里需要注意的是所有统计得到的回归方程在自变量很大或者很小的时候容易发生失真，所以回归这种模型只能用来预测和自变量统计的区间比较接近的自变量对应的函数值。

源码：《白话大数据与机器学习》