関数型プログラミング言語 F# (fsharp)の数値計算ライブラリ Math.NET Numericsを使う(7) OptimizationのLevenberg-Marquardt法（LevenbergMarquardtMinimizer）で非線形最小二乗法（回帰）でNISTの例題Rat43を計算する。

今回はこの例題。

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

コードはこんな感じで。

open MathNet.Numerics.LinearAlgebra

open MathNet.Numerics.Optimization

open Plotly.NET

let xdata = vector [ 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]

let ydata = vector [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]

let InitialValue = vector [100.0; 10.0; 1.0; 1.0]

let OptimizeFunction (p : Vector<float>) (x : Vector<float>) =

    vector [for i in 0..x.Count-1 -> p[0] / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3])]

let Gradient (p : Vector<float>) (x : Vector<float>) =

    let m =

        matrix [

            [for i in 0..x.Count-1 -> 1.0 / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3])]

            [for i in 0..x.Count-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))]

            [for i in 0..x.Count-1 -> 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))]

            [for i in 0..x.Count-1 ->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] ))]

        ]

    m.Transpose()

let obj = ObjectiveFunction.NonlinearModel(OptimizeFunction, Gradient, xdata, ydata)

let solver = LevenbergMarquardtMinimizer()

let result = solver.FindMinimum(obj, InitialValue)

let x = vector [ for i in 0..100 -> 1.0 + (float)i * 14.0/100.0 ]

let y = OptimizeFunction result.MinimizingPoint  x

[Chart.Point(xdata, ydata, ShowLegend = true, Name = "Data");

Chart.Line(x,y, ShowLegend = true, Name = "Nonlinear Fitting")]

|> Chart.combine |> Chart.withXAxisStyle("x") |> Chart.withYAxisStyle("y")

|> Chart.withTitle("Levenberg-Marquard Method") |> Chart.withLayoutStyle(Width = 800, Height = 500)

結果はこちら。