Report 3

最小二乘法

第十题希望对给定数据做最小二乘，将数据画成散点图后发现其近似一个二次函数，因此分别用 2, 3, 4 次的多项式函数拟合。

为了处理多项式，计算多项式在某点的值， well print 一个多项式，求多项式的导数与极值，我写了一个 class 抽象。以下是带有单元测试的代码。注意到其中用了 newton_solve 这个函数，这是我的第一篇报告中写到的求根的函数，这里用来处理导数的零点。

import numpy as np
from functools import reduce
from newton import newton_solve


class Polynomial(object):

    def __init__(self, coeficient):
        self.coef = coeficient

    def eval(self, x):
        f = lambda y, z: y * x + z
        res = reduce(f, self.coef[::-1])
        return res

    def derivative(self):
        new_coef = np.array([i * j for i, j in enumerate(self.coef) if i != 0])
        return Polynomial(new_coef)

    def derivative_at_x0(self, x0):
        return self.derivative().eval(x0)

    def print_poly(self):
        st = ''
        l = len(self.coef)
        # print(self.coef)
        for i, co in enumerate(self.coef[::-1][:-1]):
            st += '{}x^{}+'.format(np.array([co])[0], l-i-1)
        st += str(self.coef[0])        
        print(st)

    def extremum(self, start_point=5):
        # deri = lambda x0:self.derivative_at_x0(x0)
        x, _ = newton_solve(self.derivative_at_x0, start_point, self.derivative().derivative_at_x0)
        print('find a extremum:x={}, f(x)={}'.format(x, self.eval(x)))
        return x, self.eval(x)


def main():
    # represent polynomial: 3x^2+2x-1
    coef = [-1, 2, 3]
    temp = Polynomial(coef)
    # supposed to be 5
    # temp.eval(3)
    temp.print_poly()
    temp.derivative().print_poly()
    # print(temp.eval(5), temp.derivative_at_x0(5))
    temp.extremum()

if __name__ == '__main__':
    main()

利用 Polynomial 的抽象，我写了一个求 least square 的 class ，对于特定的数据与拟合次数，可以给出一个 Polynomial 对象。代码如下：

import numpy as np
import matplotlib.pyplot as plt
import itertools
from Polynomial import Polynomial


class LS(object):

    def __init__(self, x, y):
        self.x = x
        self.y = y

    def regression(self, degree):
        coef = self._get_coef(degree)
        return Polynomial(coef)

    def _get_coef(self, degree):
        mat = np.array([self.x**n for n in range(degree)])
        A = np.dot(mat, mat.T)
        b = np.dot(mat, np.reshape(self.y, (-1, 1)))
        coef = np.linalg.solve(A, b)
        return coef


def main():
    x_data = np.array([1, 3, 4, 5, 6, 7, 8, 9, 10])
    y_data = np.array([10, 5, 4, 2, 1, 1, 2, 3, 4])
    deg_li = [3, 4, 5]
    fig = plt.figure(figsize=(8, 8))
    ax = fig.add_subplot(111)
    marker = itertools.cycle(('+', '.', 'o', '*'))
    color = itertools.cycle(('r', 'b', 'g', 'p'))
    for i in deg_li:
        tp = LS(x_data, y_data)
        P = tp.regression(i)
        P.print_poly()
        x = np.linspace(0, 10)
        y = P.eval(x)
        ax.plot(x, y, linestyle='', marker=next(marker), markersize=10, label=r'$degree={}$'.format(i-1), color='r')
        x, y = P.extremum()
        # ax.plot(x, y, color=next(color), markersize=10)
    plt.legend(fontsize=15, frameon=False)
    plt.show()

if __name__ == '__main__':
    main()

同时单元测试里将 degree=2, 3, 4 次的情况画了出来，结果如图。

对应的最低点分别为：

find a extremum:x=[ 6.73711635], f(x)=[ 1.31496943]
find a extremum:x=[ 6.80500347], f(x)=[ 1.25338469]
find a extremum:x=[ 6.51218311], f(x)=[ 1.12282884]

Compound Simpson Method

没什么好说的，写了一个函数。

代码：

import numpy as np

def CompoundSimpson(f, a, b, eps=1e-4):
    arr = np.arange(a+eps, b, eps)
    arr2 = arr-0.5*eps
    res = eps/6*(f(a)+f(b) + 4*np.sum(f(arr2))+2*np.sum(f(arr)))
    return res
def main():
    f = lambda x:np.sqrt(1+np.cos(x)**2)
    a = 0
    b = 48
    print(CompoundSimpson(f, a, b))
if __name__ == '__main__':
    main()

输出：

58.4703899973