python 最小二乘法库_利用python搞机器学习——最小二乘法 | 学步园

原创：lhdgriver

博客：machine_love_lsearning (最好访问我的博客，因为csdn的代码高亮和公式支持不好)

邮箱：machinelovelearning@gmail.com

想了半天为啥叫最小二乘法，Least Square Method也不带这么翻译的吧。。

最小二乘法其实可以看成是一个最优化问题(好像什么问题都可以看成是最优化问题，只要你能形式化。)

一个简单的最小二乘法的例子就是Polynomial Curve Fitting(多项式曲线拟合)，就是用多项式来拟合一条曲线。转化成数学语言就是：我们有一条需要拟合的曲线$$t(x)$$，现在我们观察到了这条曲线上的n个点：

$$[(x_1, t(x_1)), (x_2, t(x_2)...(x_n, t(x_n))]$$

接下来我们要用一条多项式曲线

$$y(x,w) = w_0+w_{1}x+w_{2}x^2+...+w_{M}x^M$$

来拟合$$t(x)$$

最小二乘就是要最小化误差函数：

$$E(w)=\sum_{i=1}^{n}(y(x_i,w)-t(x_i))^2$$

python提供的库里面恰好就有做最小二乘的库，我们的$$t(x)$$选为$$sin(2 \pi x)$$，并加上一个正态分布的小噪音干扰。然后用多项式分布去拟合。

# -*- coding: utf-8 -*-

#作者:lhdgriver

#博客:lhdgriver.cn

#邮箱:machinelovelearning@gmail.com

import numpy as np #惯例

import scipy as sp #惯例

from scipy.optimize import leastsq #这里就是我们要使用的最小二乘的函数

import pylab as pl

m = 9 #多项式的次数

def real_func(x):

return np.sin(2*np.pi*x) #sin(2 pi x)

def fake_func(p, x):

f = np.poly1d(p) #多项式分布的函数

return f(x)

#残差函数

def residuals(p, y, x):

return y - fake_func(p, x)

#随机选了9个点，作为x

x = np.linspace(0, 1, 9)

#画图的时候需要的“连续”的很多个点

x_show = np.linspace(0, 1, 1000)

y0 = real_func(x)

#加入正态分布噪音后的y

y1 = [np.random.normal(0, 0.1) + y for y in y0]

#先随机产生一组多项式分布的参数

p0 = np.random.randn(m)

plsq = leastsq(residuals, p0, args=(y1, x))

print 'Fitting Parameters ：', plsq[0] #输出拟合参数

pl.plot(x_show, real_func(x_show), label='real')

pl.plot(x_show, fake_func(plsq[0], x_show), label='fitted curve')

pl.plot(x, y1, 'bo', label='with noise')

pl.legend()

pl.show()

输出的拟合参数为：

Fitting Parameters ： [ 6.33068152e+03 -2.44670137e+04 3.84443373e+04 -3.16226823e+04 1.46220505e+04 -3.77117680e+03 4.80788520e+02 -1.68548519e+01

-1.83073005e-01]

输出的图像为：

然后很明显，绿色的线过拟合了。有心的读者应该早就发现了，以上例子是Pattern Reconnition and Machine Learning上的例子。看过这本书的人都应该还记得，接下来我们应该是要加上penalty term来控制过拟合的情况。误差函数变为了：

$$E(w)=\sum_{i=1}^{n}(y(x_i,w)-t(x_i))^2+\lambda\|w\|^2$$

这个时候我们只需要改变下残差函数就行了：

#残差函数

def residuals(p, y, x):

ret = y - fake_func(p, x)

ret = np.append(ret, np.sqrt(regularization)*p) #将lambda^(1/2)p加在了返回的array的后面

return ret

设置好正则化系数后，图右：

很明显，我们的正则化约束起作用了。