最適化計算１ 直線の当てはめ 実験編

ということで，下記で理論をまとめたので直線当てはめの実験です．
daily-tech.hatenablog.com

ステップとしては７つあって，下記のステップに沿って実装してます．x

０．直線の係数定義
１．プロット用のデータを生成する．
２．フィッティングする誤差混入データを生成する．
３．生成したデータに対して，正規分布ベースの誤差をのせます．
４．最小二乗法を適用するマトリクス M を計算します．
５．固有値問題を解いて，パラメータベクトルを求めます．
６．最小二乗法で求まった直線（推定直線）を使ってデータを生成します．
７．推定直線のデータをプロットします．

ということで，結果は下記のようになりました．青が正しい直線，緑が推定直線です．当然ながら，誤差の分散を大きくしていくと青と緑の直線が合わなくなってきます．

※一点実装上の注意として！
２日くらいはまったのですが，Ax + By + C = 0 で求める係数のスケールが大きく違う (C >> A や C >> B, C << A や C << B) というような状況になると，丸め誤差の影響が顕著になってきてうまく推定直線が求まりません．これを回避するために，M行列を作成するときにスケーリングを行っています．

誤差分散 0.5

誤差分散 10

下記，実験に使った全コードです．

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
import numpy as np

def generateVecFromLine(a, b, c, minVal, maxVal, period):
  nCols = (maxVal - minVal) / period + 1
  
  if abs(a) > abs(b):
    yVec = np.linspace(minVal, maxVal, nCols)
    yVec = np.reshape(yVec, (yVec.shape[0], 1))    
    xVec = - (b / a) * yVec - (c / a)
  else:
    xVec = np.linspace(minVal, maxVal, nCols)
    xVec = np.reshape(xVec, (xVec.shape[0], 1))
    yVec = - (a / b) * xVec - (c / b) 
  
  return np.matrix(np.concatenate((xVec, yVec), axis=1))

def addGaussError(data, avg, stdv):
  dataWError = data + np.random.normal(avg, stdv, data.shape)
  return dataWError
  
def generateMMatrix(data):
  
  xLen = max(data[:, 0])[0, 0] - min(data[:, 0])[0, 0]
  yLen = max(data[:, 1])[0, 0] - min(data[:, 1])[0, 0]
  factor = max(xLen, yLen)
  
  a = np.matrix(factor * np.ones((data.shape[0], 1)))
  P = np.concatenate((data, a), axis = 1)
  return P.T * P, factor
  
def calcLineParamsFromEigenVector(M):
  lamdas, v = np.linalg.eigh(M)
  print(v[:, np.argmin(lamdas)])
  return v[:, np.argmin(lamdas)]
                      
if __name__ == "__main__":
  print("Line Fitting Sample")

  # 0. 直線の係数定義
  A = -0.01
  B = 0.025
  C = 1
  
  # 1. 直線上のデータを生成する．（真値）
  dataOrg = generateVecFromLine(A, B, C, -100, 100, 0.01)

  # 2. 誤差混入用データを生成する．（真値）
  dataSampled = generateVecFromLine(A, B, C, -50, 50, 10.0)  
  fig, ax = plt.subplots(ncols = 1, figsize=(10, 10))
  plt.xlim(-100, 100)
  plt.ylim(-100, 100)
  ax.plot(dataOrg[:, 0], dataOrg[:, 1])
  
  # 3. 真値に対して，誤差を付与する．
  dataWError = addGaussError(dataSampled, 0.0, 5.0)
  ax.scatter(dataWError[:, 0], dataWError[:, 1], s = 2, edgecolors="red")
  
  # 4. 行列 M を作成する．
  M, factor = generateMMatrix(dataWError)
  
  # 5. 固有値問題を解く
  vEig = calcLineParamsFromEigenVector(M)
   
  # 6. 最小二乗法で求まった直線からデータ生成
  dataEstEig = generateVecFromLine(vEig[0, 0], vEig[1, 0], vEig[2, 0] * factor, -100, 100, 0.01)
  
  # 7. 推定直線をプロット
  ax.plot(dataEstEig[:, 0], dataEstEig[:, 1])
  
  fig.show()

  print("vEig : ")
  print(vEig)

  print("M : ")
  print(M)
  
  print("Factor : ")
  print(factor)