c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合

科学史上众星云集,璨若星河。这些牛人基本上都是天才,但也不乏无名之辈凭借匪夷所思、骇世惊俗的猜想而跻身于巨星之列。比如,门捷列夫,整了一张留空的元素周期表,引得全世界的化学家去做填空题。还有一位德国的中学老师,名唤约翰·提丢斯(Johann Daniel Titius)的,在1766年写下了这么一个数列:

(0+4)/10 = 0.4(3+4)/10 = 0.7(6+4)/10 = 1.0(12+4)/10 = 1.6(24+4)/10 = 2.8(48+4)/10 = 5.2(96+4)/10 = 10.0(192+4)/10 = 19.6(384+4)/10 = 38.8...

当时,人们已知太阳系有六大行星,即水星、金星、地球、火星、木星、土星。如果以日地距离(约1.5亿公里)为一个天文单位,则六大行星距离太阳的距离,恰好接近提丢斯的这个数列,并且留下了无限的遐想!这个数列被称为提丢斯——波得定则。

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

1801年新年的晚上,意大利神父朱塞普·皮亚齐还在聚精会神地观察着星空。突然,他从望远镜里发现了一颗非常小的星星,正好在提丢斯——波得定则中2.8的位置上。这颗行星在几天的观测期内不断变动位置。当皮亚齐再想进一步观察这颗小行星时,他却病倒了。等到他恢复健康,再想寻找这颗小行星时,它却不知所踪。皮亚齐没有放弃这一机会,他认为这可能就是人们一直没有发现的那颗行星。 天文学家对皮亚齐的这一发现持有不同的看法。有人认为皮亚齐是正确的,也有人认为这可能是一颗彗星,关于这颗行星的说法不在少数,天文学界议论纷纷。 此时,又一位大神出现了,他就是数学天才高斯。根据皮亚齐的观测资料,高斯以其卓越的数学才能,只用了一个小时就算出了这颗神秘小行星的轨道形状,并指出它将于何时出现在哪一片天空里。1801年12月31日夜,德国天文爱好者奥伯斯,在高斯预言的时间里,用望远镜对准了这片天空。不出所料,这颗神秘小行星再一次奇迹般地出现了! 从约翰·提丢斯,到神父朱塞普·皮亚齐,再到数学王子高斯、天文爱好者奥伯斯,众多牛人联手发现的这个天体,被命名为谷神星。谷神星是太阳系中最小的、也是唯一位于小行星带的矮行星。谷神星曾经被当作一个标准:凡是比它大的行星,可以视为和地球平级的行星,否则视为小行星。比如冥王星,比谷神星大,被封为太阳系的老九,直到2006年,因为被分类到矮行星,才退出了九大行星的行列。 那么,高斯究竟是如何计算谷神星运行轨道并作出预测的呢?其实就是使用最小二乘法拟合轨迹线。高斯使用的最小二乘法发表于1809年他的著作《天体运动论》中。 让我们冒充高斯,来一次cosplay,体会一下做牛人的感觉吧。

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

假定上图是皮亚齐从零时到十时的谷神星观测记录,每一个红点对应的横坐标表示观测时间,对应的纵坐标表示谷神星的位置。如果有一条曲线(非常接近地)经过了上图的11个观测点,只要找到这条曲线的方程:

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

我们就可以随意预测11时、12时的谷神星位置了。当年高斯就是这么想的。可是,如何找到这条曲线的方程呢?这个问题难不住高斯。他用一个k次多项式 g(x)来近似 f(x):

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

通过选择合适的系数 a
0,a
1,a
2,a
3,…,a
k ,令误差:

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

为最小,就可以认为g(x)就是我们要寻找的曲线方程了。高斯的最小二乘法,就是根据观测数据和给定的多项式次数k来寻找误差最小的一组多项式系数 a0,a1,a2,a3,…,ak,有了这组系数,就可以得到 g(x),进而计算出 g(11)、g(12) 的值 ,这样就可以预测谷神星在11时、12时的位置了。

我们不是高斯,不会使用最小二乘法,幸好 numpy 给我们提供了强大的工具,使得我们可以继续冒充高斯。

>>> import numpy as np>>> _x = np.linspace(0, 10, 11)>>> _y = np.array([-0.3, -0.5, -0.2, -0.3, 0, 0.4, 0.2, -0.3, 0.2, 0.5, 0.4])>>> np.polyfit(_x, _y, 3)array([ 0.00066045, -0.01072261,  0.12684538, -0.43146853])

这里,_x 是皮亚齐的观测时间序列,_y 是皮亚齐观测到的谷神星位置序列,我们选择多项式次数为3。执行np.polyfit(xs, ys, 3),就是用最小二乘法找到三次多项式的4个系数,使得这个三次多项式和观测数据的误差最小。这个三次多项式写出来是这样的:

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

用这个函数去验证观测数据:

>>> g = np.poly1d(np.polyfit(_x, _y, 3))>>> g(_x)array([-0.43146853, -0.31468531, -0.21538462, -0.12960373, -0.05337995,        0.01724942,  0.08624709,  0.15757576,  0.23519814,  0.32307692,        0.42517483])>>> loss = np.sum(np.square(g(_x)-_y))>>> loss0.4857342657342658>>> import matplotlib.pyplot as plt>>> plt.plot(_x, _y, c='r', ls='', marker='o')>>> plt.plot(_x, g(_x), c='g', ls=':')>>> plt.show()

用三次多项式g(x)代替f(x),最小误差是0.4857。观测数据和g(x)的近似数据画在一起,效果如下

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

很显然,这个拟合效果是无法找到谷神星的。别着急,我们还可以尝试更高次幂的多项式,看看效果如何。下面的代码,同时画出了4至9次多项式的拟合结果和误差。

import numpy as npimport matplotlib.pyplot as pltplt.rcParams['font.sans-serif'] = ['FangSong']plt.rcParams['axes.unicode_minus'] = False_x = np.linspace(0, 10, 11)_y = np.array([-0.3, -0.5, -0.2, -0.3, 0, 0.4, 0.2, -0.3, 0.2, 0.5, 0.4])g3 = np.poly1d(np.polyfit(_x, _y, 3))g4 = np.poly1d(np.polyfit(_x, _y, 4))g5 = np.poly1d(np.polyfit(_x, _y, 5))g6 = np.poly1d(np.polyfit(_x, _y, 6))g7 = np.poly1d(np.polyfit(_x, _y, 7))g8 = np.poly1d(np.polyfit(_x, _y, 8))g9 = np.poly1d(np.polyfit(_x, _y, 9))loss3 = np.sum(np.square(g3(_x)-_y))loss4 = np.sum(np.square(g4(_x)-_y))loss5 = np.sum(np.square(g5(_x)-_y))loss6 = np.sum(np.square(g6(_x)-_y))loss7 = np.sum(np.square(g7(_x)-_y))loss8 = np.sum(np.square(g8(_x)-_y))loss9 = np.sum(np.square(g9(_x)-_y))plt.plot(_x, _y, c='r', ls='', marker='o')plt.plot(_x, g3(_x), label=u'三次多项式,误差%0.4f'%loss3)plt.plot(_x, g4(_x), label=u'四次多项式,误差%0.4f'%loss4)plt.plot(_x, g5(_x), label=u'五次多项式,误差%0.4f'%loss5)plt.plot(_x, g6(_x), label=u'六次多项式,误差%0.4f'%loss6)plt.plot(_x, g7(_x), label=u'七次多项式,误差%0.4f'%loss7)plt.plot(_x, g8(_x), label=u'八次多项式,误差%0.4f'%loss8)plt.plot(_x, g9(_x), label=u'九次多项式,误差%0.4f'%loss9)plt.legend()plt.show()

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

可以看到9次多项式拟合的误差,已经小到了0.0010,拟合曲线近乎精准地经过了所有的观测点。这个效果足可令人满意。下图只保留了9次多项式的拟合结果,看起来更清晰。

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

Well, cosplay is over. 真正的天体运行轨道不会这么变态,天体位置也不是一个单值就可以表示出来的。上面的例子,纯粹就是多项式拟合。有时候,待拟合数据是下图所示的样子,没关系,上面的拟合方法仍然是有效的。

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

拟合代码如下:

import numpy as npimport matplotlib.pyplot as pltplt.rcParams['font.sans-serif'] = ['FangSong']plt.rcParams['axes.unicode_minus'] = Falsexs = np.linspace(-1, 1, 11)ys = np.array([-0.3, -0.5, -0.2, -0.3, 0, 0.4, 0.2, -0.3, 0.2, 0.5, 0.4])    xm = np.linspace(-1, 1, 201)ym = ((xm**2-1)**3 + 0.5)*np.sin(2*xm) + np.random.random(201)/10 - 0.1fs4 = np.poly1d(np.polyfit(xs, ys, 4))fs5 = np.poly1d(np.polyfit(xs, ys, 5))fs6 = np.poly1d(np.polyfit(xs, ys, 6))fs7 = np.poly1d(np.polyfit(xs, ys, 7))fs8 = np.poly1d(np.polyfit(xs, ys, 8))fs9 = np.poly1d(np.polyfit(xs, ys, 9))fs10 = np.poly1d(np.polyfit(xs, ys, 10))fm2 = np.poly1d(np.polyfit(xm, ym, 2))fm3 = np.poly1d(np.polyfit(xm, ym, 3))fm4 = np.poly1d(np.polyfit(xm, ym, 4))fm5 = np.poly1d(np.polyfit(xm, ym, 5))fm6 = np.poly1d(np.polyfit(xm, ym, 6))fm7 = np.poly1d(np.polyfit(xm, ym, 7))plt.subplot(211)plt.plot(xs, ys, c='r', ls='', marker='o')plt.plot(xs, fs4(xs), label=u'四次多项式')plt.plot(xs, fs5(xs), label=u'五次多项式')plt.plot(xs, fs6(xs), label=u'六次多项式')plt.plot(xs, fs7(xs), label=u'七次多项式')plt.plot(xs, fs8(xs), label=u'八次多项式')plt.plot(xs, fs9(xs), label=u'九次多项式')plt.legend()plt.subplot(212)plt.plot(xm, ym, c='g', ls='', marker='.')plt.plot(xm, fm2(xm), label=u'二次多项式')plt.plot(xm, fm3(xm), label=u'三次多项式')plt.plot(xm, fm4(xm), label=u'四次多项式')plt.plot(xm, fm5(xm), label=u'五次多项式')plt.plot(xm, fm6(xm), label=u'六次多项式')plt.plot(xm, fm7(xm), label=u'七次多项式')plt.legend()plt.show()

拟合效果如下图所示:

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

后记

近期有很多朋友通过私信咨询有关python学习问题。为便于交流,我在CSDN的app上创建了一个小组,名为“python作业辅导小组”,面向python初学者,为大家提供咨询服务、辅导python作业。欢迎有兴趣的同学扫码加入。

《c++ 三次多项式拟合_从寻找谷神星的过程,谈最小二乘法实现多项式拟合》

    原文作者:weixin_39951396
    原文地址: https://blog.csdn.net/weixin_39951396/article/details/111049476
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞