Published on

Python で行列演算で最小二乗法，指数近似

Authors

前回 Python の Scipy を使って最小二乗法によるフィッティングを行いました．ただ最終的にベイズ線形回帰も使ってみたいと思っていたのと，中で何をしてるか追いたかったので，行列を使って最小二乗法を解くことにしました．scipy.optimize.leastsq は非線形の最小二乗法を求める用かな？

前回と同様のデータを使って，指数近似します．

y = a * exp(b * x)
log(y) = log(a) + bx

基底関数 phi を 1 + x としているから，計算した重み w は log をとったものになっているので最終的に exp をとって戻しています．

基底関数は重みの線形結合と見なすから，おそらく上記のような変形が必要．対数近似と累乗近似についても重みの線形結合に直さないといけないと思う．

下記のスクリプトを実行すると下記の重みがでて前回の結果と同じになっていることが確認できます.

0.00580943068574 4.84443568163

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

import numpy

def least_squares_regression(X, t, phi):

    PHI = numpy.array([phi(x) for x in X])
    w = numpy.linalg.solve(numpy.dot(PHI.T, PHI), numpy.dot(PHI.T, t))

    return w

if __name__ == "__main__":

    X = numpy.array([6.559112404564946264e-01,
                     6.013845740818111185e-01,
                     4.449591514877473397e-01,
                     3.557250387126167923e-01,
                     3.798882550532960423e-01,
                     3.206955701106445344e-01,
                     2.600880460776140990e-01,
                     2.245379618606005157e-01])

    t = numpy.array([1.397354195522357567e-01,
                     1.001406990711011247e-01,
                     5.173231204524778720e-02,
                     3.445520251689743879e-02,
                     3.801366557283047953e-02,
                     2.856782588754304408e-02,
                     2.036328213585812327e-02,
                     1.566228252276009869e-02])

    phi = lambda x: [1, x]

    w = least_squares_regression(X, numpy.log(t), phi)

    print numpy.exp(w[0]), w[1]

参考:

Discuss on Twitter • View on GitHub