在python中利用最小二乘拟合二次抛物线函数的方法

1、最小二乘也可以拟合二次函数

我们都知道用最小二乘拟合线性函数没有问题，那么能不能拟合二次函数甚至更高次的函数呢？答案当然是可以的。下面我们就来试试用最小二乘来拟合抛物线形状的的图像。

对于二次函数来说，一般形状为 f(x) = a*x*x+b*x+c，其中a,b,c为三个我们需要求解的参数。为了确定a、b、c，我们需要根据给定的样本，然后通过调整这些参数，知道最后找出一组参数a、b、c，使这些所有的样本点距离f(x)的距离平方和最小。用什么方法来调整这些参数呢？最常见的自然就是我们的梯度下降喽。

spicy库中有名为leastsq的方法，只需要输入一系列样本点，给出待求函数的基本形状，就可以针对上述问题求解了。

2、抛物线拟合源码

#!/usr/bin/env python
# coding:utf-8

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

# 待拟合的数据
X = np.array([1,2,3,4,5,6])
Y=np.array([9.1,18.3,32,47,69.5,94.8])

# 二次函数的标准形式
def func(params, x):
 a, b, c = params
 return a * x * x + b * x + c

# 误差函数，即拟合曲线所求的值与实际值的差
def error(params, x, y):
 return func(params, x) - y

# 对参数求解
def slovePara():
 p0 = [10, 10, 10]

 Para = leastsq(error, p0, args=(X, Y))
 return Para

# 输出最后的结果
def solution():
 Para = slovePara()
 a, b, c = Para[0]
 print "a=",a," b=",b," c=",c
 print "cost:" + str(Para[1])
 print "求解的曲线是:"
 print("y="+str(round(a,2))+"x*x+"+str(round(b,2))+"x+"+str(c))

 plt.figure(figsize=(8,6))
 plt.scatter(X, Y, color="green", label="sample data", linewidth=2)

 # 画拟合直线
 x=np.linspace(0,12,100) ##在0-15直接画100个连续点
 y=a*x*x+b*x+c ##函数式
 plt.plot(x,y,color="red",label="solution line",linewidth=2)
 plt.legend() #绘制图例
 plt.show()

solution()

上面的代码中，稍微注意的是如下几点：

1.func是待拟合的曲线的形状。本例中为二次函数的标准形式。

2.error为误差函数。很多同学会问不应该是最小平方和吗？为什么不是func(params, x) - y * func(params, x) - y？原因是名为lasts的方法中帮我们做了。看一下sklearn中源码的注释就知道什么情况了：

Minimize the sum of squares of a set of equations.
 x = arg min(sum(func(y)**2,axis=0))
   y

二次方的操作在源码中帮我们实现了。

3.p0里放的是a、b、c的初始值，这个值可以随意指定。往后随着迭代次数增加，a、b、c将会不断变化，使得error函数的值越来越小。

4.leastsq的返回值是一个tuple，它里面有两个元素，第一个元素是a、b、c的求解结果，第二个则为cost function的大小！

3.程序的最终结果与拟合曲线

程序最终的输出结果：

a= 2.06607141425 b= 2.5975001036 c= 4.68999985496
cost:1
求解的曲线是:
y=2.07x*x+2.6x+4.68999985496

最终的拟合曲线：

4、模拟其他曲线

leastsq函数除了可以模拟线性函数二次函数等多项式，还适用于任何波形的模拟。

比如方波：

def square_wave(x,p):
 a, b, c, T = p
 y = np.where(np.mod(x-b,T)<T/2, 1+c/a, 0)
 y = np.where(np.mod(x-b,T)>T/2, -1+c/a, y)
 return a*y

比如高斯分布：

def gaussian_wave(x,p):
 a, b, c, d= p
 return a*np.exp(-(x-b)**2/(2*c**2))+d

只要将上面代码中的func换成对应的函数即可！

以上这篇在python中利用最小二乘拟合二次抛物线函数的方法就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持我们。

赞 (0)

相关推荐

  • Python基于最小二乘法实现曲线拟合示例

    Python基于最小二乘法实现曲线拟合示例

    本文实例讲述了Python基于最小二乘法实现曲线拟合.分享给大家供大家参考,具体如下: 这里不手动实现最小二乘,调用scipy库中实现好的相关优化函数. 考虑如下的含有4个参数的函数式: 构造数据 import numpy as np from scipy import optimize import matplotlib.pyplot as plt def logistic4(x, A, B, C, D): return (A-D)/(1+(x/C)**B)+D def residuals(p

  • python中matplotlib实现最小二乘法拟合的过程详解

    python中matplotlib实现最小二乘法拟合的过程详解

    前言 最小二乘法Least Square Method,做为分类回归算法的基础,有着悠久的历史(由马里·勒让德于1806年提出).它通过最小化误差的平方和寻找数据的最佳函数匹配.利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小.最小二乘法还可用于曲线拟合.其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达. 下面这篇文章主要跟大家介绍了关于python中matplotlib实现最小二乘法拟合的相关内容,下面话不多说,来一起看看详细的介绍:

  • 对python实现二维函数高次拟合的示例详解

    在参加"数据挖掘"比赛中遇到了关于函数高次拟合的问题,然后就整理了一下源码,以便后期的学习与改进. 在本次"数据挖掘"比赛中感觉收获最大的还是对于神经网络的认识,在接近一周的时间里,研究了进40种神经网络模型,虽然在持续一周的挖掘比赛把自己折磨的惨不忍睹,但是收获颇丰.现在想想也挺欣慰自己在这段时间里接受新知识的能力.关于神经网络方面的理解会在后续博文中补充(刚提交完论文,还没来得及整理),先分享一下高次拟合方面的知识. # coding=utf-8 import

  • Python中实现最小二乘法思路及实现代码

    Python中实现最小二乘法思路及实现代码

    之所以说"使用"而不是"实现",是因为python的相关类库已经帮我们实现了具体算法,而我们只要学会使用就可以了.随着对技术的逐渐掌握及积累,当类库中的算法已经无法满足自身需求的时候,我们也可以尝试通过自己的方式实现各种算法. 言归正传,什么是"最小二乘法"呢? 定义:最小二乘法(又称最小平方法)是一种数学优化技术,它通过最小化误差的平方和寻找数据的最佳函数匹配. 作用:利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误

  • Python线性拟合实现函数与用法示例

    Python线性拟合实现函数与用法示例

    本文实例讲述了Python线性拟合实现函数与用法.分享给大家供大家参考,具体如下: 1. 参考别人写的: #-*- coding:utf-8 -*- import math import matplotlib.pyplot as plt def linefit(x , y): N = float(len(x)) sx,sy,sxx,syy,sxy=0,0,0,0,0 for i in range(0,int(N)): sx += x[i] sy += y[i] sxx += x[i]*x[i]

  • 在python中利用最小二乘拟合二次抛物线函数的方法

    1.最小二乘也可以拟合二次函数 我们都知道用最小二乘拟合线性函数没有问题,那么能不能拟合二次函数甚至更高次的函数呢?答案当然是可以的.下面我们就来试试用最小二乘来拟合抛物线形状的的图像. 对于二次函数来说,一般形状为 f(x) = a*x*x+b*x+c,其中a,b,c为三个我们需要求解的参数.为了确定a.b.c,我们需要根据给定的样本,然后通过调整这些参数,知道最后找出一组参数a.b.c,使这些所有的样本点距离f(x)的距离平方和最小.用什么方法来调整这些参数呢?最常见的自然就是我们的梯度下降

  • 在python中利用opencv简单做图片比对的方法

    下面代码中利用了两种比对的方法,一 对图片矩阵(m x m)求解特征值,通过比较特征值是否在一定的范围内,判断图片是否相同.二 对图片矩阵(m x m)中1求和,通过比较sum和来比较图片. # -*- coding: utf-8 -*- import cv2 as cv import numpy as np import os file_dir_a='C:\Users\wt\Desktop\data\image1\\' file_dir_b='C:\Users\wt\Desktop\data\

  • 在python中利用GDAL对tif文件进行读写的方法

    利用GDAL库对tif影像进行读取 示例代码默认波段为[B.G.R.NIR的顺序,且为四个波段] import gdal def readTif(fileName): dataset = gdal.Open(fileName) if dataset == None: print(fileName+"文件无法打开") return im_width = dataset.RasterXSize #栅格矩阵的列数 im_height = dataset.RasterYSize #栅格矩阵的行

  • 在python中利用numpy求解多项式以及多项式拟合的方法

    构建一个二阶多项式:x^2 - 4x + 3 多项式求解 >>> p = np.poly1d([1,-4,3]) #二阶多项式系数 >>> p(0) #自变量为0时多项式的值 3 >>> p.roots #多项式的根 array([3., 1.]) >>> p(p.roots) #多项式根处的值 array([0., 0.]) >>> p.order #多项式的阶数 2 >>> p.coeffs #

  • python中利用zfill方法自动给数字前面补0

    python中利用zfill方法自动给数字前面补0

    python中有一个zfill方法用来给字符串前面补0,非常有用 view sourceprint? n = "123" s = n.zfill(5) assert s == "00123" zfill()也可以给负数补0 n = "-123" s = n.zfill(5) assert s == "-0123" 对于纯数字,我们也可以通过格式化的方式来补0 n = 123 s = "%05d" % n a

  • python中利用h5py模块读取h5文件中的主键方法

    如下所示: import h5py import numpy as np #HDF5的写入: imgData = np.zeros((2,4)) f = h5py.File('HDF5_FILE.h5','w') #创建一个h5文件,文件指针是f f['data'] = imgData #将数据写入文件的主键data下面 f['labels'] = np.array([1,2,3,4,5]) #将数据写入文件的主键labels下面 f.close() #关闭文件 #HDF5的读取: f = h5

  • python中利用numpy.array()实现俩个数值列表的对应相加方法

    小编想把用python将列表[1,1,1,1,1,1,1,1,1,1] 和 列表 [2,2,2,2,2,2,2,2,2,2]对应相加成[3,3,3,3,3,3,3,3,3,3]. 代码如下: import numpy a = numpy.array([1,1,1,1,1,1,1,1,1,1]) b = numpy.array([2,2,2,2,2,2,2,2,2,2]) c = a + b print(type(c)) print(list(c)) 输出结果为: <class 'numpy.nd

  • 在python中利用KNN实现对iris进行分类的方法

    如下所示: from sklearn.datasets import load_iris iris = load_iris() print iris.data.shape from sklearn.cross_validation import train_test_split X_train, X_test, y_train, y_test = train_test_split(iris.data, iris.target, test_size = 0.25, random_state = 3

  • python中利用matplotlib读取灰度图的例子

    代码为: import matplotlib.pyplot as plt #用于显示图片 import matplotlib.image as mpimg # mpimg 用于读取图片 picdir = 'C:\\Users\\wyt\\Desktop\\test\\ship\\012400.jpg' img = mpimg.imread(picdir) plt.imshow(img) plt.title('Original train image') plt.show() 显示结果如下,是以热

  • 在python中利用try..except来代替if..else的用法

    在有些情况下,利用try-except来捕捉异常可以起到代替if-else的作用. 比如在判断一个链表是否存在环的leetcode题目中,初始代码是这样的 # Definition for singly-linked list. # class ListNode(object): # def __init__(self, x): # self.val = x # self.next = None class Solution(object): def hasCycle(self, head):

随机推荐