专栏首页机器学习与统计学Lagrange、Newton、分段插值法及Python实现

Lagrange、Newton、分段插值法及Python实现

数据分析中,经常需要根据已知的函数点进行数据、模型的处理和分析,而通常情况下现有的数据是极少的,不足以支撑分析的进行,这里就需要使用差值法模拟新的数值来满足需求。

插值法又称“内插法”,是利用函数f(x)在某区间中已知的若干点的函数值,作出适当的特定函数,在区间的其他点上用这特定函数的值作为函数f (x)的近似值。

常用的插值方法有Lagrange插值、Newton插值、分段插值、Hermite插值、样条插值等等。这里我们就介绍一下最常用到的Lagrange、Newton、分段插值法及Python实现。

1、拉格朗日插值法 Lagrange插值基本思想是将待求的n次多项式插值函数pn(x)改写成另一种表示方式,再利用插值条件确定其中的待定函数,从而求出插值多项式。它是n次多项式插值,成功地用构造插值基函数的方法解决了求n次多项式插值函数问题。

一般地,若已知

在互不相同 n+1 个点

处的函数值

( 即该函数过

这n+1个点),则可以考虑构造一个过这n+1 个点的、次数不超过n的多项式

,使其满足:

要估计任一点ξ,ξ≠xi,i=0,1,2,...,n,则可以用Pn(ξ)的值作为准确值f(ξ)的近似值。称式(*)为插值条件(准则),含xi(i=0,1,...,n)的最小区间[a,b],其中a=min{x0,x1,...,xn},b=max{x0,x1,...,xn}。

设集合

是关于点

的角标的集合,

,作n个多项式

。对于任意

,都有

使得

是n-1次多项式,且满足

并且

。最后可得

形如上式的插值多项式

称为拉格朗日(Lagrange)插值多项式。

#coding=utf-8
from matplotlib import pyplot as plt

def Lg(data,testdata):
    predict=0
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    if testdata in data_x:
        #print "testdata is already known"
        return data_y[data_x.index(testdata)]
    for i in range(len(data_x)):
        af=1
        for j in range(len(data_x)):
            if j!=i:
                af*=(1.0*(testdata-data_x[j])/(data_x[i]-data_x[j]))
        predict+=data_y[i]*af
    return predict

def plot(data,nums):
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    Area=[min(data_x),max(data_x)]
    X=[Area[0]+1.0*i*(Area[1]-Area[0])/nums for i in range(nums)]
    X[len(X)-1]=Area[1]
    Y=[Lg(data,x) for x in X]
    plt.plot(X,Y,label='result')
    for i in range(len(data_x)):
        plt.plot(data_x[i],data_y[i],'ro',label="point")
    plt.savefig('Lg.jpg')
    plt.show()

data=[[0,0],[1,2],[2,3],[3,8],[4,2]]
print Lg(data,1.5)
plot(data,100)

2、牛顿插值 Newton插值基本思想是将待求的n次插值多项式Pn(x)改写为具有承袭性的形式,然后利用插值条件⑴确定Pn(x)的待定系数,以求出所要的插值函数。

牛顿差值引入了差商的概念,使其在差值节点增加时便于计算。

设函数

,已知其n+1个插值节点为

,我们定义:

的零阶差商为

在点

的一阶差商为

在点

的二阶插商为

一般的,

在点

的k 阶差商为

可将k阶差商表示为函数值

的组合:

经过分别变形,依次代入,可得牛顿差值公式:

可记为:

其中,

为牛顿差值公式的余项或截断误差,当n趋于无穷大时为零。

#coding=utf-8
#coding=utf-8
from matplotlib import pyplot as plt

def calF(data):
    #差商计算  n个数据 0-(n-1)阶个差商 n个数据
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    F= [1 for i in range(len(data))]   
    FM=[]
    for i in range(len(data)):
        FME=[]
        if i==0:
            FME=data_y
        else:
            for j in range(len(FM[len(FM)-1])-1):
                delta=data_x[i+j]-data_x[j]
                value=1.0*(FM[len(FM)-1][j+1]-FM[len(FM)-1][j])/delta
                FME.append(value)
        FM.append(FME)
    F=[fme[0] for fme in FM]
    print FM
    return F

def NT(data,testdata,F):
    #差商之类的计算
    predict=0
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    if testdata in data_x:
        return data_y[data_x.index(testdata)]
    else:
        for i in range(len(data_x)):
            Eq=1
            if i!=0:
                for j in range(i):
                    Eq=Eq*(testdata-data_x[j])
            predict+=(F[i]*Eq)
        return predict

def plot(data,nums):
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    Area=[min(data_x),max(data_x)]
    X=[Area[0]+1.0*i*(Area[1]-Area[0])/nums for i in range(nums)]
    X[len(X)-1]=Area[1]
    F=calF(data)
    Y=[NT(data,x,F) for x in X]
    plt.plot(X,Y,label='result')
    for i in range(len(data_x)):
        plt.plot(data_x[i],data_y[i],'ro',label="point")
    plt.savefig('Newton.jpg')
    plt.show()

data=[[0,0],[1,2],[2,3],[3,8],[4,2]]

plot(data,100)

3、分段线性插值 对每一个分段区间(xi,xi+1)分别进行插值,将被插值函数f(x)的插值节点由小到大排序,然后每对相邻的两个节点为端点的区间上用m次多项式去近似f(x)。

将每两个相邻的节点用直线连起来,如此形成的一条折线就是分段线性插值函数。计算x点的插值时,只用到x左右的两个节点,计算量与节点个数n(初始值x0,y0的长度,n=length(x0))无关,而拉格朗日插值与n值有关。分段线性插值中n越大,分段越多,插值误差越小。

#coding=utf-8

from matplotlib import pyplot as plt
def DivideLine(data,testdata):
    #找到最邻近的
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    if testdata in data_x:
        return data_y[data_x.index(testdata)]
    else:
        index=0
        for j in range(len(data_x)):
            if data_x[j]<testdata and  data_x[j+1]>testdata:
                index=j
                break
        predict=1.0*(testdata-data_x[j])*(data_y[j+1]-data_y[j])/(data_x[j+1]-data_x[j])+data_y[j]
        return predict

def plot(data,nums):
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    Area=[min(data_x),max(data_x)]
    X=[Area[0]+1.0*i*(Area[1]-Area[0])/nums for i in range(nums)]
    X[len(X)-1]=Area[1]
    Y=[DivideLine(data,x) for x in X]
    plt.plot(X,Y,label='result')
    for i in range(len(data_x)):
        plt.plot(data_x[i],data_y[i],'ro',label="point")
    plt.savefig('DivLine.jpg')
    plt.show()

data=[[0,0],[1,2],[2,3],[3,8],[4,2]]
print DivideLine(data,1.5)
plot(data,100)

参考文献: http://www.cnblogs.com/duye/p/8671820.html https://www.cnblogs.com/ytxwzqin/p/9539659.html 代码引用:BUAA-XX@CSDN https://blog.csdn.net/sinat_33829806/article/details/78387843

本文分享自微信公众号 - 机器学习与统计学(tjxj666),作者:统计学家

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-04-28

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 数据科学家易犯的十大编码错误,你中招了吗?

    我是一名高级数据科学家,在 Stackoverflow 的 python 编码中排前 1%,而且还与众多(初级)数据科学家一起工作。下文列出了我常见到的 10 ...

    统计学家
  • 【Python基础系列】常见的数据预处理方法(附代码)

    本文简单介绍python中一些常见的数据预处理,包括数据加载、缺失值处理、异常值处理、描述性变量转换为数值型、训练集测试集划分、数据规范化。

    统计学家
  • 【数据分析 R语言实战】学习笔记 第三章 数据预处理 (上)

    数据是分析的核心,在做数据分析之前,首先要对数据进行一定的处理。数据预处理指当录入或读取数据后,对数据进行必要的清理,包括查错纠错、异常观察值和无效样本的处理、...

    统计学家
  • 为什么你应该学习 Python 的生成器?

    写过一段时间代码的同学,应该对这一句话深有体会:程序的时间利用率和空间利用率往往是矛盾的,可以用时间换空间,可以用空间换时间,但很难同时提高一个程序的时间利用率...

    青南
  • 框架篇-Vue面试题1-为什么 vue 组件中的 data 是函数而不是对象

    当一个组件被定义,data必须声明为返回一个初始数据对象的函数,因为组件可能被用来创建多个实例

    itclanCoder
  • Python实现删除某列中含有空值的行的示例代码

    砸漏
  • C++ string实现

    作为C++从业者,我相信都会被考察过实现简单的string类,包括构造、析构、拷贝构造以及赋值拷贝等,因为这能够很好的考察面试者的C++基本功。借看《剑指off...

    evenleo
  • Python中时间格式数据的处理

    1、时间转换 时间转换是指字符型的时间格式数据,转换成为时间型数据的过程。 一般从csv导入过来的文件,时间都保存为字符型格式的,需要转换。 时间转换函数: d...

    Erin
  • R分类算法-神经网络算法

    神经网络(Artifical Neural Network) 神经网络(人工神经网络),是一种模仿生物网络(动物的中枢神经系统,特别是大脑)的结构和功能的数学模...

    Erin
  • R分类算法-Logistic回归算法

    逻辑回归 Logistic Regression 所谓LR,就是一个被Logistic方程归一化后的线性回归,可以将非线性的问题转化为线性问题。 优点: ...

    Erin

扫码关注云+社区

领取腾讯云代金券