【Hello Python World】Week 12：Matplotlib：Plot gracefully

画出函数的图像

分析

这道题是对matplotlib中的一些简单方法的应用，我尝试了：
* 修改函数曲线的颜色，线型，线宽
* 调整坐标轴的限制
* 设置坐标轴、图例信息
* 获得函数的在某点的极值(使用argmax)并标注出来
* 设置垂直标线，以及在图上加入一个特殊标注
Python在画图上功能也比较强大，个人感觉用Python画图和`MATLAB`画图的区别和
用`Word`与用`Latex`写文档一样，入门较难，熟悉以后会比另一种更方便一些。

代码

# coding: utf-8
from matplotlib.pyplot import *
from numpy import *

y = []
for x in linspace(0, 2, 2000):
    y.append((sin(x-2))**2*exp(-x**2))

plot(linspace(0, 2, 2000), y, 
    label = '$y = sin^2 (x-2) e^{-x^2}$', 
    linewidth=2.0, linestyle="-", color='m')
    #颜色是洋红色，线型为实线，线宽2.0

#对坐标轴的限制
xlim(0, 2)
ylim(0) #从0截起

title('Exercise 11.1 Plot a function') #标题
xlabel('$x$', color = 'b', fontsize = 16) #设置坐标轴信息
ylabel('$y = sin^2 (x-2) e^{-x^2}$',color = 'b', fontsize = 16)
legend(loc = 'upper right') #显示图例在右上角
grid() #显示网格

t = linspace(0, 2, 2000)[argmax(y)] #获得最高这一点的x轴坐标
plot([t,t],[0,(sin(t-2))**2*exp(-t**2)],
    color = 'r', linewidth = 3, linestyle = ':') #画出标记的红色虚线
scatter(t,(sin(t-2))**2*exp(-t**2), 50, color ='r') #标出最高点

annotate(r'$x_{max}=%f$' %t, color = 'g',
         xy=(t, 0), xycoords='data', #xy表示的元组是针对的地址
         xytext=(20, 40), textcoords='offset points', fontsize=16,
         arrowprops=dict(color = 'g', arrowstyle="->", connectionstyle="arc3,rad=.4"))

show()

输出结果如下

参数估计

分析

关于公式$\boldsymbol{\hat{b}}=\arg\min_{\boldsymbol{b}}\left \| X\boldsymbol{b}-\boldsymbol{y} \right \|_2$，相信不少人会和我一开始想的一样，想要用Python中的argmin方法来解决。但是这个地方并不适合使用这些统计的方法，因为没有一个可供筛选的数据集啊！什么意思呢？想想一下使用argmin或者min的时候，有一个参数就是要有可迭代的数据集(比如列表)，答案已经在这个数据集里了，这样解更像是一道选择题，而在这里我们是不知道答案是什么的，甚至连范围都不知道(你也许会说：b不是已经知道了吗？可是那是拿来验证的啊，不能当做已知条件！)，解决起来更像是一道填空题，答案需要自己推导(想用穷举法？星辰大海，从何开始？)。
求解的方法是最小二乘法，这是多元线性回归里最常用的方法之一，定义可见wiki，我们直接使用结论：b估计解为$\hat {\boldsymbol{b}}=(X^{T}X)^{-1}X^{T}{\boldsymbol {y}}$,，具体的推导过程可见这篇博客，推导过程详尽完整。
当然，Python中也有提供相应的方法可以直接求解最小二乘问题，方法是scipy.optimize.leastsq，方法原型

scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, 
                        col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08,
                        gtol=0.0, maxfev=0, epsfcn=None, factor=100, 
                        diag=None)

一般我们只要指定前三个参数就可以了，在本问题中使用的形式是leastsq(error, b0, (X,y))，error是自定义的估计误差的函数，第二个元素是寻找b的迭代起始点b0，返回值是一个二元组，第一个元素是即是b，第二个元素用不到。完整的用法可见官方文档。

代码

# coding: utf-8
from numpy import *
from scipy.optimize import *
from scipy.linalg import *
from matplotlib.pyplot import *

X = random.randint(-5, 5, size=(20,10))
b = np.random.randint(-5, 5, size=10)
z = np.random.normal(size=20)
y = np.dot(X,b) + z

#方法1：公式法
e_b = dot(dot(inv(dot(X.T, X)), X.T), y)

#方法2：使用现成方法
def error(b_e, X, y):
    #return norm(dot(X, b_e)-y, 2) 错误！只需要做差就行了，不需要把完整的范数写出来
    return dot(X, b_e)-y

b0 = random.random(10) #需要指定一个迭代的初始点
ret = leastsq(error, b0, (X,y)) 
e_b = ret[0] #取二元组的第一个参数就是解了
#--------------------------------------------------

#开始画图
title('Exercise 11.2 Estimate the solution of a linear system') #标题

xlim(0, 11)
ylim(-5.5, 5.5)
#将自己的符号与刻度值形成一一映射，替换掉原有的数值刻度，使用列表解析的方法生成带不同下标的符号
xticks(range(1, 11), ['$b_' + str(i) + '$' for i in range(0, 10)]) 
xlabel('Estimator of $b_k$', color = 'b', fontsize = 16) #设置坐标轴信息
ylabel('Value',color = 'b', fontsize = 16)

scatter(range(1, 11), b , 30, color ='limegreen', label = 'True value of $b$')
scatter(range(1, 11), e_b , 25, color ='r', label = 'Estimator of $b$', marker = '*')

hlines(0,0,11 ,colors = 'b', linestyles = ':') #画出y=0那条虚线
legend(loc = 'upper right') #显示图例在右上角
grid()
show()

结果如图所示

直方图和密度估计

分析

写到这里，可以发现这些题都是和统计或机器学习有关系的，这些知识在数学课上也没有涉及到，所以关于方法也比较陌生，这个题考察的是核密度估计。为了完成这道题，我们需要：1）画概率分布直方图；2) 用高斯核密度估计法求密度估计(数学原理见wiki或博客)
关于核密度估计，scipy.stats提供了相应的方法gaussian_kde用来求核，详情见文档，核密度估计的意思是根据生成的数据可以估计出一个核(类似于曲线拟合里面根据数据拟合出的一个函数)，用例kernel = gaussian_kde(sample)。

代码

# coding: utf-8
from numpy import *
from scipy.stats import *
from matplotlib.pyplot import *

observe_data = random.beta(1.5, 1.2, 10000) #随机生成10000个贝塔分布的观测值
#根据这些观测值，我们可以用高斯核密度估计来获得这个分布的核
kernel = gaussian_kde(observe_data) 

title('Exercise 11.3 Histogram and Gauss KDE with $\\beta$-Distribution') #标题
xlabel('x', color = 'b', fontsize = 16) #设置坐标轴信息
ylabel('Observer value',color = 'b', fontsize = 16)

hist(observe_data, bins = 25, color = 'g', alpha = 0.3, 
    density = True, label = 'Probability density curve') #画直方图
plot(arange(0, 1, 0.0001), kernel.pdf(arange(0, 1, 0.0001)), 
    color = 'r', label = 'Probability density curve') #画密度曲线
grid()
legend()
show()

结果如下