Wikipedia見ながら書いたらちょっとバグったけど持ち直したラグランジュ補間

基本は前回と同じ
https://note.com/alchan/n/nc14a3c453144

Java(processing)の場合

        //始点ー終点間を補間
       ArrayList<PVector> LagrangeInterpolation(ArrayList<PVector> vectors, float t)
       {
           PVector sv = vectors.get(0);
           PVector ev = vectors.get(vectors.size() - 1);

           float total_x = ev.x - sv.x;
           //float dx = total_x*t;

           ArrayList<PVector> ret = new ArrayList<PVector>();

           ret.add(vectors.get(0).copy());
           for (float i = t; i < 1; i += t)
           {
               float x = (float)(sv.x + (total_x * i));
               float y = (float)(Lx(vectors, x));
               ret.add(new PVector(x, y));
           }
           ret.add(vectors.get(vectors.size() - 1).copy());
           return ret;
       }
       //ラグランジュ補間(一回分のf(x))
       double Lx(ArrayList<PVector> vectors, double x)
       {
           double ret = 0;
           for (int i = 0; i < vectors.size(); i++)
           {
               ret += vectors.get(i).y * ljx(vectors, i, x);
           }
           return ret;
       }
       //ラグランジュ基底関数
       double ljx(ArrayList<PVector> vectors, int j, double x)
       {
           double ret = 1;
           
           //jまででなく、vectorsの最後まで回す
           //for (int i = 0; i < j; i++)            
           for (int i = 0; i < vectors.size(); i++)
           {
               if (i == j) { continue; }
               double numerator = x - vectors.get(i).x;
               double denominator = vectors.get(j).x - vectors.get(i).x;
               ret *= (numerator / denominator);
           }
           return ret;
       }

Pythonの場合

Google Colabで確認

#ラグランジュ基底関
def ljx(vx, j, x):
  ret = 1
  for i in range(len(vx)):
  #for i, item in enumerate(vx):
     if i == j:
       continue
     numerator = x - vx[i]
     denominator = vx[j] - vx[i]
     ret *= (numerator / denominator)
  return ret
#ラグランジュ補間(一回分のf(x))
def Lx(vx, vy, x):
  ret = 0
  for i in range(len(vx)):
    ret += vy[i] * ljx(vx, i, x)
  return ret
#始点ー終点間を補間
def LagrangeInterpolation(vx, vy, t):
  total_x = vx[-1]-vx[0]  
  rx=[]
  ry=[]
  for i in np.arange(0,1,t):
    x = vx[0]+total_x*i
    y = Lx(vx,vy,x)
    rx.append(x)
    ry.append(y)
  rx.append(vx[-1])
  ry.append(vy[-1])
  return (rx,ry)

確認用コード

import numpy as np
import matplotlib.pyplot as plt

#直線のパラメータ
a = 1
b = 0
#正規分布ノイズの振れ幅
deviation = 1

#ラグランジュ基底関
def ljx(vx, j, x):
  ret = 1
  for i in range(len(vx)):
  #for i, item in enumerate(vx):
     if i == j:
       continue
     numerator = x - vx[i]
     denominator = vx[j] - vx[i]
     ret *= (numerator / denominator)
  return ret
#ラグランジュ補間(一回分のf(x))
def Lx(vx, vy, x):
  ret = 0
  for i in range(len(vx)):
    ret += vy[i] * ljx(vx, i, x)
  return ret
#始点ー終点間を補間
def LagrangeInterpolation(vx, vy, t):
  total_x = vx[-1]-vx[0]  
  rx=[]
  ry=[]
  for i in np.arange(0,1,t):
    x = vx[0]+total_x*i
    y = Lx(vx,vy,x)
    rx.append(x)
    ry.append(y)
  rx.append(vx[-1])
  ry.append(vy[-1])
  return (rx,ry)

#正規分布のノイズをのせた直線
def f(x):
  r = np.random.randn(x.size)#正規分布
  y = a*x+(deviation*r)+b
  return y

#上限と下限を設定できる乱数列
def get_random_sequence(count, min, max):
  seq = (max-min)*np.random.rand(count)+min
  return seq

x0 = get_random_sequence(20, 1, 5)
x0.sort()
y0 = f(x0)

v = LagrangeInterpolation(x0,y0,0.01)
vx = v[0]
vy = v[1]

plt.plot(vx, vy, color="k") 
plt.scatter(x0, y0)

plt.title('lagrange')
plt.show()

a = 1
b = 0
をパラメータとする直線に

偏差deviation = 1
で正規分布に従う乱数をのせたもの

その散布図

これをラグランジュ補間で接続することを考える。
結果

y値の振れが大き過ぎて散布図を極度に圧縮する結果を得たと思われる。