计算pi的近似值公式法_Python金融大数据分析：金融学中最常用的数学技术之一逼近法...

首先，是通常的导入工作：

In [1]: import numpy as np
        from pylab import plt, mpl
In [2]: plt.style.use('seaborn')
        mpl.rcParams['font.family'] = 'serif'
        %matplotlib inline

本节使用的主函数示例如下，由一个三角函数项和一个线性项组成：

In [3]: def f(x):
            return np.sin(x) + 0.5 * x

重点是在给定区间内通过回归和插值求取该函数的近似值。首先，生成该函数的图形，以便更好地观察逼近法的效果。我们感兴趣的区间是[−2π，2π]。图11-1显示了该函数在np.linspace()函数定义的固定区间上的图像。create_plot()是一个助手函数，可以创建本章多次要使用的同类图表：

图11-1　示例函数图表

In [4]: def create_plot(x, y, styles, labels, axlabels):
            plt.figure(figsize=(10, 6))
            for i in range(len(x)):
                plt.plot(x[i], y[i], styles[i], label=labels[i])
                plt.xlabel(axlabels[0])
                plt.ylabel(axlabels[1])
            plt.legend(loc=0)
In [5]: x = np.linspace(-2 * np.pi, 2 * np.pi, 50) ❶
In [6]: create_plot([x], [f(x)], ['b'], ['f(x)'], ['x', 'f(x)'])

❶ 用于绘图和计算的x值。

11.1.1　回归

回归是相当高效的函数近似值计算工具。它不仅适用于求取一维函数的近似值，在更高维度上也很有效。得出回归结果所需要的数值化方法很容易实现，执行也很快速。本质上，回归的任务是在给定一组所谓“基函数”bd，d∈{1，…，D}的情况下，根据公式11-1找出最优参数

，…

，其中对于i∈{1,…I}观察点，yi ≡ f(xi)。xi可以视为自变量观测值，yi可视为因变量观测值(从函数或者统计的意义上说)。

公式11-1.　最小化回归问题

1．作为基函数的单项式

最简单的情况是以单项式作为基函数——也就是说，b1=1，b2=x，b3=x2，b4=x3，…在这种情况下，NumPy有可以确定最优参数(np.polyfit())和通过一组输入值求取近似值(np.polyval())的内建函数。

表11-1列出了np.polyfit()函数的参数。在np.polyfit()返回的最优回归相关系数ρ基础上，np.polyval(ρ,x)返回x坐标的回归值。

典型向量化风格的np.polyfit()和np.polyval()线性回归(deg=1)的应用方式如下。由于回归估算值保存在ry数组中，所以我们可以如图11-2那样比较回归结果和原始函数。当然，线性回归无法处理示例函数的sin部分：

In [7]: res = np.polyfit(x, f(x), deg=1, full=True) ❶
In [8]: res ❷
Out[8]: (array([ 4.28841952e-01, -1.31499950e-16]),
         array([21.03238686]),
         2,
         array([1., 1.]),
         1.1102230246251565e-14)
In [9]: ry = np.polyval(res[0], x) ❸
In [10]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
                     ['f(x)', 'regression'], ['x', 'f(x)'])

❶ 线性回归步骤。

❷ 完整的结果：回归参数、残差、有效秩、奇异值和相对条件数。

❸ 使用回归参数求值。

图11-2　线性回归

为了处理示例函数的sin部分，必须使用更高次的单项式。下一个回归试图使用5次单项式作为基函数。果不其然，回归结果(如图11-3所示)看上去更接近原始函数。但是，它还远称不上完美：

In [11]: reg = np.polyfit(x, f(x), deg=5)
         ry = np.polyval(reg, x)
In [12]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
                     ['f(x)', 'regression'], ['x', 'f(x)'])

图11-3　使用最高5次的单项式进行回归

最后一次尝试使用7次的单项式作为基函数来计算示例函数的近似值。这次的结果如图11-4所示，相当有说服力：

图11-4　7次单项式回归

In [13]: reg = np.polyfit(x, f(x), 7)
         ry = np.polyval(reg, x)
In [14]: np.allclose(f(x), ry)❶
Out[14]: False
In [15]: np.mean((f(x) - ry) ** 2) ❷
Out[15]: 0.0017769134759517689
In [16]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
                     ['f(x)', 'regression'], ['x', 'f(x)'])

❶ 检查函数和回归值是否相同(至少接近)。

❷ 根据函数值计算回归值均方差(MSE)。

2．单独的基函数

一般来说，当您选择更好的基函数组时，可以得到更好的回归结果，例如利用对函数的认识进行近似值计算。在这种情况下，单独的基函数必须通过一个矩阵方法定义(也就是使用NumPy的ndarray对象)。首先，例子中的多项式最高为3次(图11-5)。本例的核心函数是np.linalg.lstsq()：

In [17]: matrix = np.zeros((3 + 1, len(x))) ❶
         matrix[3, :] = x ** 3 ❷
         matrix[2, :] = x ** 2 ❷
         matrix[1, :] = x ❷
         matrix[0, :] = 1 ❷
In [18]: reg = np.linalg.lstsq(matrix.T, f(x), rcond=None)[0] ❸
In [19]: reg.round(4) ❹
Out[19]: array([ 0. , 0.5628, -0. , -0.0054])
In [20]: ry = np.dot(reg, matrix) ❺
In [21]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
                     ['f(x)', 'regression'], ['x', 'f(x)'])

❶ 基函数值(矩阵)所用的ndarray对象。

❷ 从常数到三次基函数值。

❸ 回归步骤。

❹ 最优回归参数。

❺ 函数值的回归估算。

图11-5　有单独基函数的回归

根据前面单项式的经验，图11-5中的结果并不真的如预期那么好。使用更通用的方法可以让我们利用对示例函数的认识。我们知道函数中有一个sin部分。因此，在基函数中包含一个正弦函数是有意义的。简单起见，我们替换最高次的单项式。现在的拟合很完美，如图11-6所示：

In [22]: matrix[3, :] = np.sin(x) ❶
In [23]: reg = np.linalg.lstsq(matrix.T, f(x), rcond=None)[0]
In [24]: reg.round(4) ❷
Out[24]: array([0. , 0.5, 0. , 1. ])
In [25]: ry = np.dot(reg, matrix)
In [26]: np.allclose(f(x), ry) ❸
Out[26]: True
In [27]: np.mean((f(x) - ry) ** 2) ❸
Out[27]: 3.404735992885531e-31
In [28]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
                     ['f(x)', 'regression'], ['x', 'f(x)'])

❶ 新的基函数利用关于示例函数的知识。

❷ 最优回归参数恢复原始参数。

❸ 现在，回归产生了完美的拟合。

图11-6　使用正弦基函数的回归

3．有噪声的数据

回归对于有噪声的数据同样能够很好的处理，这种数据来自于模拟或者(不完善的)测量。为了阐述这个要点，我们生成同样具有噪声的自变量观测值和因变量观测值。图11-7表明，回归结果比有噪声的数据点更接近原始函数。在某种意义上，回归在一定程度上平均了噪声：

图11-7　使用有噪声数据的回归

In [29]: xn = np.linspace(-2 * np.pi, 2 * np.pi, 50) ❶
         xn = xn + 0.15 * np.random.standard_normal(len(xn)) ❷
         yn = f(xn) + 0.25 * np.random.standard_normal(len(xn)) ❸
In [30]: reg = np.polyfit(xn, yn, 7)
         ry = np.polyval(reg, xn)
In [31]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
                     ['f(x)', 'regression'], ['x', 'f(x)'])

❶ 新的x确定值。

❷ 在x值中引入噪声。

❸ 在y值中引入噪声。

4．未排序数据

回归的另一个重要特点是，它可以无缝地处理未排序数据。前面的例子都依赖于经过排序的x数据，情况并不总是这样的。为了说明这一点，我们随机生成自变量数据点。在这种情况下，仅靠从视觉上检查原始数据很难识别出任何结构：

In [32]: xu = np.random.rand(50) * 4 * np.pi - 2 * np.pi ❶
         yu = f(xu)
In [33]: print(xu[:10].round(2)) ❶
         print(yu[:10].round(2)) ❶
         [-4.17 -0.11 -1.91 2.33 3.34 -0.96 5.81 4.92 -4.56 -5.42]
         [-1.23 -0.17 -1.9 1.89 1.47 -1.29 2.45 1.48 -1.29 -1.95]
In [34]: reg = np.polyfit(xu, yu, 5)
         ry = np.polyval(reg, xu)
In [35]: create_plot([xu, xu], [yu, ry], ['b.', 'ro'],
                     ['f(x)', 'regression'], ['x', 'f(x)'])

❶ 随机化x值。

和有噪声数据一样，回归方法不关心观测点的顺序。这在研究公式11-1所示的最小化问题的结构时很明显。从图11-8中显示的结果来看也很明显。

图11-8　使用未排序数据的回归

5．多维

最小二乘回归方法的另一个优点是，不需要太多的修改就可以用于多维的情况。接下来以fm()函数为例进行讲解：

In [36]: def fm(p):
             x, y = p
             return np.sin(x) + 0.25 * x + np.sqrt(y) + 0.05 * y ** 2

为了正确地可视化这个函数，我们需要自变量数据点的网格(在两个维度上)。图11-9根据以x、y、z表示的自变量和因变量二维数据点网格，显示了fm()函数的形状：

In [37]: x = np.linspace(0, 10, 20)
         y = np.linspace(0, 10, 20)
         X, Y = np.meshgrid(x, y) ❶
In [38]: Z = fm((X, Y))
         x = X.flatten() ❷
         y = Y.flatten() ❷
In [39]: from mpl_toolkits.mplot3d import Axes3D ❸
In [40]: fig = plt.figure(figsize=(10, 6))
         ax = fig.gca(projection='3d')
         surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2,
                                cmap='coolwarm', linewidth=0.5,
                                antialiased=True)
         ax.set_xlabel('x')
         ax.set_ylabel('y')
         ax.set_zlabel('f(x, y)')
         fig.colorbar(surf, shrink=0.5, aspect=5)

❶ 从一维ndarray对象生成二维ndarray对象(网格)。

❷ 从二维ndarray对象得到一维ndarray对象。

❸ 必要时从matplotlib导入3D绘图功能。

图11-9　使用两个参数的函数

为了获得好的回归结果，我们将应用基本函数集，包括fm()函数、np.sin()和np.sqrt()函数。图11-10直观展示了完美的回归结果：

In [41]: matrix = np.zeros((len(x), 6 + 1))
         matrix[:, 6] = np.sqrt(y) ❶
         matrix[:, 5] = np.sin(x) ❷
         matrix[:, 4] = y ** 2
         matrix[:, 3] = x ** 2
         matrix[:, 2] = y
         matrix[:, 1] = x
         matrix[:, 0] = 1
In [42]: reg = np.linalg.lstsq(matrix, fm((x, y)), rcond=None)[0]
In [43]: RZ = np.dot(matrix, reg).reshape((20, 20)) ❸
In [44]: fig = plt.figure(figsize=(10, 6))
         ax = fig.gca(projection='3d')
         surf1 = ax.plot_surface(X, Y, Z, rstride=2, cstride=2,
                     cmap=mpl.cm.coolwarm, linewidth=0.5,
                     antialiased=True) ❹
         surf2 = ax.plot_wireframe(X, Y, RZ, rstride=2, cstride=2,
                                   label='regression') ❺
         ax.set_xlabel('x')
         ax.set_ylabel('y')
         ax.set_zlabel('f(x, y)')
         ax.legend()
         fig.colorbar(surf, shrink=0.5, aspect=5)

❶ 用于y参数的np.sqrt()函数。

❷ 用于x参数的np.sin()函数。

❸ 将回归结果转化为网格结构。

❹ 绘制原始函数曲面。

❺ 绘制回归曲面。

图11-10　双参数函数的回归曲面

回归

最小二乘回归方法有多种应用领域，包括简单的函数逼近和基于有噪声或者未排序数据的函数逼近。这些方法可以应用于一维问题，也可以应用于多维问题。由于这种方法的基础数学理论，所以它在一维问题和多维问题上的应用总是“几乎相同”。

11.1.2　插值

与回归相比，插值(例如，3次样条插值)在数学上更为复杂。它还被限制在低维度问题上。给定一组有序的观测点(按照x维排序)，基本的思路是在两个相邻数据点之间进行回归，这样做不仅产生的分段插值函数完全匹配数据点，而且函数在数据点上连续可微分。连续可微分性需要至少3阶插值——也就是3次样条插值。然而，这种方法一般也适用于4次或者线性样条插值。

下面的代码可实现线性样条插值，结果如图11-11所示：

In [45]: import scipy.interpolate as spi ❶
In [46]: x = np.linspace(-2 * np.pi, 2 * np.pi, 25)
In [47]: def f(x):
             return np.sin(x) + 0.5 * x
In [48]: ipo = spi.splrep(x, f(x), k=1) ❷
In [49]: iy = spi.splev(x, ipo) ❸
In [50]: np.allclose(f(x), iy) ❹
Out[50]: True
In [51]: create_plot([x, x], [f(x), iy], ['b', 'ro'],
                     ['f(x)', 'interpolation'], ['x', 'f(x)'])

❶ 从SciPy导入必要的子库。

❷ 实现线性样条插值。

❸ 得出内插值。

❹ 检查内插值是否(足够)接近函数值。

图11-11　线性样条插值(完整数据集)

如果有按照x值排序的一组数据点，那么应用本身也和使用np.polyfit()和np.polyval()函数一样简单。在本例中，它们对应的函数是sci.splrep()和sci.splev()。表11-2列出了sci.splrep()函数的主要参数。

表11-3列出sci.splev()函数的参数。

样条插值在金融学中往往用于估算未包含在原始观测点中的自变量数据点的因变量值。为此，在下个例子中选择一个更小的区间，仔细观察一次样条插入的值。图11-12说明，插值函数确实可以线性地在两个观测点之间插值。对于某些应用，这可能不够精确。此外，很明显函数在原始数据点上不是连续可微分的——这是另一个不足：

In [52]: xd = np.linspace(1.0, 3.0, 50) ❶
         iyd = spi.splev(xd, ipo)
In [53]: create_plot([xd, xd], [f(xd), iyd], ['b', 'ro'],
                     ['f(x)', 'interpolation'], ['x', 'f(x)'])

❶ 具有更多数据点的较小区间。

图11-12　线性样条插值(数据子集)

重复整个练习，这次使用3次样条，结果明显改善(见图11-13)：

In [54]: ipo = spi.splrep(x, f(x), k=3) ❶
         iyd = spi.splev(xd, ipo) ❷
In [55]: np.allclose(f(xd), iyd) ❸
Out[55]: False
In [56]: np.mean((f(xd) - iyd) ** 2) ❹
Out[56]: 1.1349319851436892e-08
In [57]: create_plot([xd, xd], [f(xd), iyd], ['b', 'ro'],
                     ['f(x)', 'interpolation'], ['x', 'f(x)'])

❶ 完整数据集上的3次样条插值。

❷ 结果应用到更小的时间间隔。

❸ 插值仍然不完美。

❹ 但好于从前。

插值

在可以应用样条插值的情况下，可以期望得到比最小二乘回归方法更好的近似结果。但是要记住，必须有排序(且“无噪声”)的数据，该方法仅限于低维度问题。样条插值的计算要求也更高，在某些用例中可能导致花费的时间比回归方法的长得多。

图11-13　三次样条插值(数据子集)

这是一本书中的内容 你知道是哪本吗？需要此书电子书的加群：1136192749