`

拉格朗日插值-python

 
阅读更多

在数据库中,有些数据是异常值或者空值,这些值在分析的时候应该特殊处理,比如最简单的忽略掉或者通过算法推测它的值。其中拉格朗日插值就是通过其他已经知道的值,对x位置缺失的值插入的算法。

 

假定我们已经知道了(x0,y0),(x1,y1),(x2,y2),(x3,y3) 四个数据点,如果有多项式L(x)刚好穿过这4个点,这样的公式可写成:




当x=x0时候,L(x0)= y0*1 + y1*0 + y2*0 + y3*0 = y0
 当x=x1时候,L(x1)= y0*0 + y1*1 + y2*0 + y3*0 = y1

...

故,L(x)刚好穿过(x0,y0),(x1,y1),(x2,y2),(x3,y3) 四个点。

只要求到函数l0,l1,l2, l3,就可以得到L(x), 求剩余点的值。

 

对于函数l0,它要满足两个两个条件:

1)当x!= x0时候,l0(x0) = 0。即x 为x1,x2,x3时候 l0(x0) = 0, 下面这个多项式即可满足该条件:

(x-x1)(x-x2)(x-x3)

2)当x=x0时, l0(x0) = 1. 可以看出在上面的多项式(x-x1)(x-x2)(x-x3)中,代入x=x0,得到一个数字k = (x0-x1)(x0-x2)(x0-x3),为了让多项式(x-x1)(x-x2)(x-x3)也满足x=x1时,它的值为1,那么该多项式除以数字k即可。

这样最终得到

l0(x) = (x-x1)(x-x2)(x-x3) / (x0-x1)(x0-x2)(x0-x3)

 

同理得到:

l1(x) = (x-x0)(x-x2)(x-x3) / (x1-x0)(x0-x2)(x0-x3)

l2(x) = (x-x0)(x-x1)(x-x3) / (x2-x0)(x2-x1)(x0-x3)

l3(x) = (x-x0)(x-x1)(x-x2) / (x3-x0)(x3-x1)(x3-x2)

 

得到了l1(x),l2(x),l3(x),代入L(x)中,最终得到了刚好穿过(x0,y0),(x1,y1),(x2,y2),(x3,y3)四个点的一个多项式。

现在可以对x=x4, x=x5...求他们的y值了。

 

下面是python的一个例子。我先生成给出了10个点,y=x*x, x = {-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10},

通过这十个点,得到一个刚好穿过这10个点的拉格朗日函数,最后通过该函数,补全了x = {-9, -7, -5...7, 9}上的值。

左边是补全之前,右边是补全了数据之后:



 

 

import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import lagrange

axis_x = [i for i in range(-10, 10, 2)]
axis_y = [i ** 2 for i in axis_x]
data = pd.Series(axis_y, index=axis_x)
plt.figure("lagrange demo")
ax1 = plt.subplot(121)
plt.sca(ax1)
plt.plot(axis_x, axis_y, linestyle=' ', marker='o', color='b')
plt.plot(axis_x, axis_y, linestyle='--', color='r')
##generate lagrange function
lag = lagrange(axis_x, axis_y)

new_x = [i for i in range(-10, 10, 1)]
new_y = []
for x in new_x:
    if x in data.index:
        new_y.append(data[x])
    else:
        new_y.append(lag(x))


ax2 = plt.subplot(122)
plt.sca(ax2)
plt.plot(new_x, new_y, linestyle=' ', marker='o', color='b')
plt.plot(new_x, new_y, linestyle='--', color='r')

plt.show()

 

 

 

  • 大小: 110.2 KB
  • 大小: 50.1 KB
分享到:
评论

相关推荐

Global site tag (gtag.js) - Google Analytics