文章目录

近来做试验报告，学分子动力学的时候处理了一些数据，已经忘得差不多的数值计算方法又捡起来一些。然而问题总是层出不穷的，现在记下来，以后再遇到的时候不至于手足无措。

关于python

窃以为用 python 做数值分析比 matlab 要好，因为相比语法简单而有简易的 matlab， python 是一门完备的语言，并且在机器学习方面正在大放异彩。再说用盗版 matlab 算出来的东西不好意思往论文上写……

Numpy

可以说 Numpy 的应用占据了 python 数据分析的一半内容，我另外开了一篇博文作为笔记。

Jupyter Notebook

python 最好的交互界面就是 Jupyter Notebook，但是这东西的配置比较麻烦。

常用魔术语句

将 numpy 和 matplotlib.pylab 的名字空间完全导入，并设置图像输出在 Notebook 中：

1
2
3

%pylab inline
%xmode Plain
%config InlineBackend.figure_format="svg"

引入 sympy 并设置默认使用 pprint ：

1 2	import sympy as sp sp.init_printing(use_latex=True)

配置 matplotlib 的字体

matplotlib 的默认字体并不支持中文，令人恼火，后来我在 segmentfault 找到了解决方案。

进入 C:/Anaconda/Lib/site-packages/matplotlib/mpl-data 修改 matplotlibrc 配置文件。搜索 font.serif ，在字体列表里加上你效用的字体名称。我改成了这样：

#font.size           : 12.0
#font.serif          : Droid Sans Fallback, New Century Schoolbook, Century Schoolbook L, 
                       Utopia, ITC Bookman, Bookman, Nimbus Roman No9 L, Times New Roman,
                       Times, Palatino, Charter, serif
font.sans-serif      : Droid Sans Fallback, Lucida Grande, Verdana, Geneva, Lucid, Arial,
                       Helvetica, Avant Garde, sans-serif

也就是使用了 Sans 字体，把首选换成了字符集全面的 Droid Sans Fallback。

之后进入 ./fonts/ttf ，删掉所有 Vera 系列的字体，并且把我们的 DroidSans.ttf 复制过来。
另外负号显示会有问题，把 axes.unicode_minus 项改成 False 即可。

Sympy 中的 mathjax 渲染

默认使用 unicode 字符画，经过之前的设置会启用 mathjax。但是默认字体看着不舒服，最后在 StackOverflow 找到了解决方案。从 mathjax 官网下载原装包，解压之后：

把 jax/output/HTML-CSS/fonts/TeX 复制到 ../notebook/static/components/MathJax/jax/output/HTML-CSS/fonts/ 目录下
把 fonts/HTML-CSS/TeX 复制到 ../notebook/static/components/MathJax/fonts/HTML-CSS/ 目录下
打开 ../notebook/static/notebook/js/main.min.js, 搜索 availableFonts，把代码修改成：

availableFonts: ["STIX-Web","TeX"],
imageFont: null;
preferredFont: "TeX",
webFont: "TeX"

不过即使修改好了还是会有一些问题，例如字号问题导致上下标一团糟。可以在查看矩阵输出的时候把 mathjax 的设置临时改一下增大字号。

神经网络库——Keras

之前研究过一段时间，经历了“自己编写\(\to\)使用 TensorFlow \(\to\)使用 Keras 作为前端”的过程。Keras 的官方教程已经有大佬翻译，在这里仅写出最常用的。

引入函数

常用函数主要在 keras.models 和 keras.layers ：

from keras.models import Sequential, load_model
from keras.layers.core import Dense, Dropout, Reshape
from keras.layers.convolutional import Conv2D
from keras.layers.pooling import MaxPooling2D

模型建立

一个最简单的序贯模型，仅需要把所用的网络层逐个写上去：

model=Sequential([                                        # 序贯模型
        Reshape([28, 28, 1], input_shape=(784,)),         # 变形层
        Conv2D(32, 5, activation='relu', padding='same'), # 2D卷积层
        MaxPooling2D(),                                   # 2D最大值池化层
        Conv2D(64, 5, activation='relu', padding='same'),  
        MaxPooling2D(), 
        Reshape([7 * 7 * 64]), 
        Dense(units=1024, activation='relu'),             # 全连接层
        Dropout(0.1),                                     # 随机断开层
        Dense(units=10, activation='softmax')
    ])

Reshape( 输出变量形态 )
Conv2D( 卷积核数目, 卷积核边长, activation='激励函数', padding='输出数据与原来的大小关系' )
MaxPooling2D()
Dense( units=神经元数目, activation='激励函数' )
Dropout( 随机断开概率 )
第一层一般加上 input_shape=输入变量形态 以方便检查变量是否合法

模型训练

# 模型编译
model.compile(optimizer='rmsprop',             # 优化器
             loss='categorical_crossentropy',  # 损失函数
             metrics=['accuracy'])             # 目标（想要最小化的）变量

# 模型训练
history=model.fit(mnist.train.images,          # 训练集.数据
                 mnist.train.labels,           # 训练集.标准输出
                 batch_size=100,               # 每个训练周期取的训练样本量
                 epochs=3)                     # 重复训练次数

模型应用

# 模型评估
model.evaluate(mnist.test.images,   # 测试集.数据
              mnist.test.labels)    # 测试集.标准输出

# 进行预测
model.predict(data)                 # 输入数据

关于Gnuplot

Gnuplot 非常方便——对于我这种渣渣机械硬盘，开启速度比 matlab 快十倍。但是作为命令行程序，没有那么直观，并且配置也比较奇怪。

Ubuntu子系统下的 Gnuplot

为了配合 LAMMPS，最好把 Gnuplot 装在 Ubuntu子系统里，但是这就有个 GUI 的问题。可以使用Windows 下的 X11 服务：

在 Windows 下安装 Xming，启动
在 Linux 下安装x11-apps
（非必须）检验 Gnuplot 的版本，安装 gnuplot-x11
使用 DISPLAY=:0.0 gnuplot 启动，如果显示终端是 wxt 的话，使用 set term x11
为了省事，可以直接在 ~/.bashrc 里加上 export DISPLAY=:0.0

结果差强人意，x11下的输出确实不敢恭维，但是仅用来看看已经够了，设置好之后再用 set term svg ，set output filename 换 svg 输出。

其中某一次，出于未知的原因（反正忘记 set term x11 了），我成功的打开了 gnuplot 的 wxt 终端（而且这个窗口的图标还是Xming），绘图质量甚至好于 Windows 下的 Gnuplot（Windows的辣鸡字体渲染……），但是之后再没成功过。

于是现在使用的方式是打开 cmd，在 bash 里用 LAMMPS 算好，再回到 cmd 画图……

输出图片的设置

如果设置输出为 png，有一些选项可以设定：

1	set term png [size x,y] [font font_name,font_size]

看起来可以随意设置输出图片的大小，获得高分辨率的图像。然而实际上真的只是调整了大小，分辨率是不变的，输出接近于一张白纸……所以还是不要用 png 输出了。

可以使用 svg，矢量图吼啊！而且可以在 Inkscape 进一步编辑，加上各种辅助线，甚至调节颜色，进一步美化。
如果要包含公式，可以考虑使用 epslatex 输出，不过也可以用 svg 输出再后期添上去。

关于算法

插值算法

一维二阶插值

等间距插值比较简单，教程也多，而通式则比较复杂。

以二阶插值为例，在各数据点之外还需要一个边界条件，一般是第一点处的导数，可以设定成 0 或者前两个点连线的斜率。那这样就可以给出插值参数的递推公式。

设点 \((x_i,y_i)\) 与 \((x_{i+1},y_{i+1})\) 之间的曲线 \(f_i(x)\) 为 \(y=a_i x^2+b_i x+c_i\) ，则 \(f_i(x)\) 满足：

\[ \left\{ \begin{split} y_i&=a_i x_i^2+b_i x_i+c_i\\ y_{i+1}&=a_i x_{i+1}^2+b_i x_{i+1}+c_i\\ 2a_i x_i+b_i&=2a_{i-1}x_i+b_{i-1}\Leftrightarrow f_{i-1}'(x_i)=f_{i}'(x_i) \end{split} \right. \]

以 \(\vec{k}=(a_i\;b_i\;c_i)^T\) 为未知量，列出矩阵方程 \(\mathbf{A}\vec{k}=\vec{b}\) ，其中：

\[ \mathbf{A}=\left( \begin{matrix} x_i^2 & x_i & 1\\ x_{i+1}^2 & x_{i+1} & 1\\ 2x_i & 1 & 0 \end{matrix} \right), \vec{b}=\left( \begin{matrix} y_i \\ y_{i+1}\\ 2a_{i-1}x_i+y_{i-1} \end{matrix} \right) \]

这样就可以在 python 中利用 numpy 很方便地写出插值函数。

二维一阶插值

高维插值在有限元中十分有用。对于三角形面元，可以比较简单的表达为行列式形式：

\[ \left( \begin{matrix} T_i \\ T_j\\ T_k \end{matrix} \right) =\left( \begin{matrix} 1&x_i&y_i\\ 1&x_j&y_j\\ 1&x_k&y_k \end{matrix} \right) \left( \begin{matrix} a_1\\ a_2\\ a_3 \end{matrix} \right) \]

很容易求出系数\((a_1\ a_2\ a_3)^T\)。

求导算法

数值积分的算法多，numpy 也有内置函数，而数值求导却不好找，于是我自己写了一个，这里以二阶导数为例。

数值求导的原理就是泰勒展开，将 \(f(x)\) 展开到二阶可得：

\[ f(x_0+t)=f(x_0)+f'(x_0)t+\frac{1}{2}f''(x_0)t^2+O(t^2) \]

可以看出，如果要以 \(O(t^2)\) 的精确度得到 \(f''(x_0)\) ，要列三个方程，因此对于三个样本点 \((x_0-t_{-1},y_{-1}),(x_0,y_0),(x_0+t_1,y_1)\)，得出三个展开式。目的是\(Ay_{-1}+By_0+Cy_1=f''(x_0)\) ，则以 \(\vec{k}=(A\;B\;C)^T\) 为未知量忽略高阶小量，列出矩阵方程 \(\mathbf{A}\vec{k}=\vec{b}\)，其中

\[ \mathbf{A}=\left( \begin{matrix} 1 & 1 & 1 \\ -t\_{-1}&0&t_1 \\ \frac{1}{2}t_{-1}^2 & 0 & \frac{1}{2}t_1^2 \end{matrix} \right), \vec{b}=\left( \begin{matrix} 0\\ 0\\ 1 \end{matrix} \right) \]

这里有一个特殊情况：当 \(t_{-1}=t_1\) 时，矩阵奇异，因此应跳过矩阵方程直接令 \(\vec{k}=(1 \ -2\ 1)^T\)。同样的，可以用 numpy 写出函数：

import numpy as np

def diff2(num):
    length=num.shape[0]
    f2=zeros([length,1])
    for i in arange(1,length-1):
        t_n,t_p=num[i,0]-num[i-1,0],num[i+1,0]-num[i,0]
        if t_n==0 or t_p==0:
            factor=matrix([[0],[0],[0]])
        elif t_n == t_p:
            factor=matrix([[1],[-2],[1]])/t_n**2
        else:
            A=matrix([[1,1,1],[t_p,0,-t_n],[t_p**2/2,0,t_n**2/2]])
            b=matrix([[0],[0],[1]])
            factor=A**-1*b
        f2[i,0]=num[i-1:i+2,1].T*factor
    f2[0,0],f2[-1,0]=f2[1,0],f2[-2,0]
    return c_[num[:,0],f2]

在实际应用之后我终于明白为什么数值求导不常用了：尽管理论上拥有 \(O(t^2)\) 的精确度，实际上每个值还是只由三个点决定。对于一般的实验信号，其中都会存在高频噪音，这种算法得出的只是在 0 附近剧烈震荡的曲线，其振幅远超曲线总体变化的二阶导数，掩盖了有用信息。因此在数值求导之前必须进行滤波，以均值滤波为佳。在我所用的例子中，滤波窗口达到约 100 时才得到了满意的结果。

统计与回归

统计与回归在实验数据处理时非常有用……但我总是用完就忘了。

总体数量特征与样本数量特征

数学期望\(\mu=E(X)\)、方差\(\sigma^2=D(X)\)都是针对随机事件的本质（总体）的，而采用样本代表总体时会有偏差。

样本均值\(\bar X=\frac 1 n \sum X_i\)，数学期望即为总体均值，即\(E(\bar X)=\mu\)，方差\(D(\bar X)=\frac{\sigma^2}{n}\)

样本方差\(S^2=\frac{1}{n-1}\sum(X_i-\bar X)^2\)，数学期望即为总体方差，即\(E(S^2)=\sigma^2\)。注意此处分母为\(n-1\)，定性的来说，这是因为算式中使用了\(X_i-\bar X\)而非\(X_i-\mu\)，\(\bar X\)本身就存在误差，因而误差增大，使得\(\frac{1}{n}\sum(X_i-\mu)^2=\frac{1}{n-1}\sum(X_i-\bar X)^2\)。

在Numpy中，有两个相关函数：

cov(X,Y)返回一个协方差矩阵，只给一个参数的话返回\(\frac{1}{n-1} \sum(X_i-\bar X)^2\)
var(X)返回一组数据的方差，即\(\frac 1 n \sum(X_i-\bar X)^2\)

拟合优度和线性相关系数

拟合优度（又称决定系数）\(R^2\)代表“模型能在多大程度上解释因变量的变化”，用于评价拟合的好坏。

以\(\bar{y}\)表示样本均值，\(\hat{y}\)表示根据回归模型预测的值，则残差平方和与方差之比越小，表明能解释的部分越多，即 \[ R^2=1-\frac{\mathrm{var}(\varepsilon)}{\mathrm{var}(Y)}=1-\frac{\sum(y_i-\hat{y})^2}{\sum(y_i-\bar{y})^2} \] 可以看出\(R^2\in(-\infty,1]\)，但一般而言是不会小于0的。

线性相关系数表明两个变量间线性相关的程度。总体相关系数\(\rho\)定义为两变量协方差与其标准差之积的比，即 \[ \rho_{X,Y}=\frac{\mathrm{cov}(X,Y)}{\sigma_X \sigma_Y} \] 那么 \(\rho\in[-1,1]\)，\(\rho=1\)表示正线性关系，\(\rho=-1\)表示负线性关系，\(\rho=0\)表示无线性关系。

实际上需要用样本估计总体，样本相关系数 \[ r=\frac{\sum(X_i-\overline{X})(Y_i-\overline{Y})}{\sqrt{\sum(X_i-\overline{X})^2}\sqrt{\sum(Y_i-\overline{Y})^2}} \] 对于简单的线性模型，\(R^2=r^2\)。

写成python代码的形式：

1
2
3

import numpy as np
R_square=lambda yi,y_:1-sum((yi-y_)**2)/(sum((yi-np.mean(yi))**2))
rho=lambda X,Y:np.cov(X,Y)[0,1]/np.sqrt(np.cov(X,Y)[0,0]*np.cov(X,Y)[1,1])

埃雅仁迪尔之星

Luke的个人博客

数据分析中的一些问题