UnityでVisual C#用の数値計算ライブラリMath.NET Numericsを使う(7) OptimizationのLevenberg-Marquardt法（LevenbergMarquardtMinimizer）で非線形最小二乗法（回帰）でNISTの例題Rat43を計算してその軌跡を描く。

さて今回は非線形最小二乗法。

こちらを移植したもの。

 Visual C# (C_sharp)の数値計算ライブラリ MathNET Numericsを使う(7)　OptimizationのLevenberg-Marquardt法（LevenbergMarquardtMinimizer）で非線形最小二乗法（回帰）でNISTの例題Rat43を計算する。

using UnityEngine;

using MathNet.Numerics;

using MathNet.Numerics.LinearAlgebra;

using MathNet.Numerics.Optimization;

using System;

public class MathNETController : MonoBehaviour

{

    public GameObject PointsPrefab;

    GameObject line;

    Vector<double> x;

    Vector<double> y;

    float scale;

    int n;

    int count;

    // Start is called once before the first execution of Update after the MonoBehaviour is created

    void Start()

    {

        var obj = ObjectiveFunction.NonlinearModel(OptimizeFunction, Gradient, Xdata, Ydata);

        var solver = new LevenbergMarquardtMinimizer();

        var result = solver.FindMinimum(obj, InitialValue1);

        x = Vector<double>.Build.DenseOfArray(Generate.LinearSpaced(1000, 1.0, 15.0));

        y = OptimizeFunction(result.MinimizingPoint, x);

        line = GameObject.Find("Line");

        scale = 50.0f;

        n = x.Count;

        count = 0;

        for (int i = 0; i < Xdata.Count; i++)

        {

            GameObject point = Instantiate(PointsPrefab);

            point.transform.position = new Vector3((float)Xdata[i], Convert.ToSingle(Ydata[i]/scale), 0f);

        }

    }

    // Update is called once per frame

    void Update()

    {

        if (count < n) {

            line.transform.position = new Vector3((float)x[count], Convert.ToSingle(y[count] / scale), 0f);

            count++;

        }

    }

    Vector<double> Xdata = Vector<double>.Build.DenseOfArray(new double[] {

            1.00,   2.00,   3.00,   4.00,   5.00,   6.00,   7.00,   8.00,   9.00,   10.00,

            11.00,  12.00,  13.00,  14.00,  15.00

        });

    Vector<double> Ydata = Vector<double>.Build.DenseOfArray(new double[] {

             16.08, 33.83,  65.80,  97.20,  191.55, 326.20, 386.87, 520.53, 590.03, 651.92,

            724.93, 699.56, 689.96, 637.56, 717.41

        });

    Vector<double> InitialValue1 = Vector<double>.Build.DenseOfArray(new double[] { 100, 10, 1, 1 });

    Vector<double> OptimizeFunction(Vector<double> p, Vector<double> x)

    {

        var y = Vector<double>.Build.Dense(x.Count);

        for (int i = 0; i < x.Count; i++)

        {

            y[i] = p[0] / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]);

        }

        return y;

    }

    private Matrix<double> Gradient(Vector<double> p, Vector<double> x)

    {

        var prime = Matrix<double>.Build.Dense(x.Count, p.Count);

        for (int i = 0; i < x.Count; i++)

        {

            prime[i, 0] = 1.0 / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]);

            prime[i, 1] = -p[0] * Math.Exp(p[1] - p[2] * x[i]) /

                (p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3] + 1.0));

            prime[i, 2] = p[0] * x[i] * Math.Exp(p[1] - p[2] * x[i]) /

                (p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3] + 1.0));

            prime[i, 3] = Math.Log(1.0 + Math.Exp(p[1] - p[2] * x[i])) * p[0] /

                (p[3] * p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]));

        }

        return prime;

    }

}

結果の動画。計算結果を玉が軌跡を描きながら動く…