Clojure(JVMで動くLISP系）でJavaの数値計算ライブラリApache Commons Mathを使う(7) OptimizationのLevenberg-Marquardt法（LevenbergMarquardtOptimizer）で非線形最小二乗法（回帰）でNISTの例題Rat43を計算する。

今回はこの例題。

Javaの数値計算ライブラリApache Commons Mathを使う(7)　OptimizationのLevenberg-Marquardt法（LevenbergMarquardtOptimizer）で非線形最小二乗法（回帰）でNISTの例題Rat43を計算する。Google Gemini 2.5 Proに助けてもらいながら…

これもChatGPTに助けてもらってすぐできた。コードはこちら。

(ns mini.NistRat43Fitter

  (:import

   [javax.swing JFrame]

   [java.awt BorderLayout]

   [org.jfree.chart ChartFactory ChartPanel]

   [org.jfree.chart.plot PlotOrientation]

   [org.jfree.chart.renderer.xy XYLineAndShapeRenderer]

   [org.jfree.data.xy XYSeries XYSeriesCollection]

   [org.apache.commons.math3.fitting.leastsquares LeastSquaresBuilder

    LevenbergMarquardtOptimizer

    MultivariateJacobianFunction]

   [org.apache.commons.math3.linear Array2DRowRealMatrix

    ArrayRealVector]

   [org.apache.commons.math3.util Pair]))

;; -----------------------------

;; モデル関数 (Rat43)

;; -----------------------------

(defn model-function [xdata]

  (reify MultivariateJacobianFunction

    (value [_ params]

      (let [b1 (.getEntry params 0)

            b2 (.getEntry params 1)

            b3 (.getEntry params 2)

            b4 (.getEntry params 3)

            n  (alength xdata)

            values   (ArrayRealVector. n)

            jacobian (Array2DRowRealMatrix. n 4)]

        (dotimes [i n]

          (let [xi (aget xdata i)

                expval (Math/exp (- b2 (* b3 xi)))

                denom (Math/pow (inc expval) (/ 1.0 b4))

                f (/ b1 denom)]

            (.setEntry values i f)

            ;; 偏微分

            (.setEntry jacobian i 0 (/ 1.0 denom))

            (.setEntry jacobian i 1

                       (- (/ (* b1 expval)

                             (* b4 (Math/pow (inc expval) (+ 1.0 (/ 1.0 b4)))))))

            (.setEntry jacobian i 2

                       (/ (* b1 xi expval)

                          (* b4 (Math/pow (inc expval) (+ 1.0 (/ 1.0 b4))))))

            (.setEntry jacobian i 3

                       (/ (* (Math/log (inc expval)) b1)

                          (* b4 b4 denom)))))

        (Pair. values jacobian)))))

;; -----------------------------

;; フィット処理

;; -----------------------------

(defn fit-data []

  (let [xdata (double-array [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])

        yobs  (double-array [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])

        start (double-array [100.0 10.0 1.0 1.0])

        problem (.. (LeastSquaresBuilder.)

                    (start start)

                    (model (model-function xdata))

                    (target yobs)

                    (lazyEvaluation false)

                    (maxEvaluations 1000)

                    (maxIterations 1000)

                    build)

        optimizer (LevenbergMarquardtOptimizer.)

        optimum (.optimize optimizer problem)

        params  (.getPoint optimum)]

    (println "NIST Rat43 (Ratkowsky3) 最小二乗フィット")

    (println "反復回数:" (.getIterations optimum))

    (println "評価回数:" (.getEvaluations optimum))

    (doseq [i (range 4)]

      (println (str "b" (inc i) ": " (.getEntry params i))))

    {:xdata xdata :yobs yobs :params params :optimum optimum}))

;; -----------------------------

;; fitted function

;; -----------------------------

(defn fitted-fn [x b1 b2 b3 b4]

  (/ b1 (Math/pow (inc (Math/exp (- b2 (* b3 x))))

                  (/ 1.0 b4))))

;; -----------------------------

;; グラフデータ作成

;; -----------------------------

(defn create-dataset [{:keys [xdata yobs params]}]

  (let [series1 (XYSeries. "Original Points")

        series2 (XYSeries. "Curve Fitting")]

    (dotimes [i (alength xdata)]

      (.add series1 (aget xdata i) (aget yobs i)))

    (dotimes [i 100]

      (let [xi (+ 1.0 (* (- 15.0 1.0) (/ (double i) 99.0)))

            yi (fitted-fn xi

                          (.getEntry params 0)

                          (.getEntry params 1)

                          (.getEntry params 2)

                          (.getEntry params 3))]

        (.add series2 xi yi)))

    (doto (XYSeriesCollection.)

      (.addSeries series1)

      (.addSeries series2))))

;; -----------------------------

;; グラフ表示

;; -----------------------------

(defn show-chart []

  (let [fit (fit-data)

        dataset (create-dataset fit)

        chart (ChartFactory/createXYLineChart

               "Levenberg-Marquardt Optimizer"

               "x" "y"

               dataset

               PlotOrientation/VERTICAL

               true false false)

        plot (.getXYPlot chart)

        renderer (XYLineAndShapeRenderer.)]

    (.setRange (.getRangeAxis plot) 0.0 800.0)

    (.setSeriesLinesVisible renderer 0 false)

    (.setSeriesShapesVisible renderer 1 false)

    (.setRenderer plot renderer)

    (doto (JFrame. "NIST Rat43 problem")

      (.setDefaultCloseOperation JFrame/EXIT_ON_CLOSE)

      (.setBounds 10 10 640 480)

      (.add (ChartPanel. chart) BorderLayout/CENTER)

      (.setVisible true))))

;; 実行エントリーポイント

(defn -main []

  (show-chart))

(-main)

結果：