１０点のデータから最小二乗法で２次関数を求めるプログラムの作り方

実現したいこと

c#
1using System;
2
3class LeastSquaresQuadratic
4{
5    static void Main()
6    {
7        // サンプルデータ（x, y）を定義
8        double[] xData = { -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 };
9        double[] yData = { 16, 9, 4, 1, 0, 1, 4, 9, 16, 25 };
10
11        // 必要な各種合計値を初期化
12        int n = xData.Length;
13        double sumX = 0, sumY = 0, sumX2 = 0, sumX3 = 0, sumX4 = 0;
14        double sumXY = 0, sumX2Y = 0;
15
16        // 各値を計算
17        for (int i = 0; i < n; i++)
18        {
19            double x = xData[i];
20            double y = yData[i];
21            sumX += x;
22            sumY += y;
23            sumX2 += x * x;
24            sumX3 += x * x * x;
25            sumX4 += x * x * x * x;
26            sumXY += x * y;
27            sumX2Y += x * x * y;
28        }
29
30        // 行列の各要素を計算
31        double[,] matrix = {
32            { sumX4, sumX3, sumX2 },
33            { sumX3, sumX2, sumX },
34            { sumX2, sumX, n }
35        };
36        double[] vector = { sumX2Y, sumXY, sumY };
37
38        // 連立方程式を解く
39        double[] coefficients = GaussianElimination(matrix, vector);
40
41        // 結果出力
42        Console.WriteLine($"y = {coefficients[0]:F3}x^2 + {coefficients[1]:F3}x + {coefficients[2]:F3}");
43    }
44
45    // ガウス消去法で連立方程式を解くメソッド
46    static double[] GaussianElimination(double[,] matrix, double[] vector)
47    {
48        int n = vector.Length;
49
50        for (int i = 0; i < n; i++)
51        {
52            // 対角要素がゼロでないようにするためのピボット選択
53            int maxRow = i;
54            for (int k = i + 1; k < n; k++)
55            {
56                if (Math.Abs(matrix[k, i]) > Math.Abs(matrix[maxRow, i]))
57                    maxRow = k;
58            }
59
60            // 行の入れ替え
61            for (int k = i; k < n; k++)
62            {
63                double tmp = matrix[maxRow, k];
64                matrix[maxRow, k] = matrix[i, k];
65                matrix[i, k] = tmp;
66            }
67            double tmp2 = vector[maxRow];
68            vector[maxRow] = vector[i];
69            vector[i] = tmp2;
70
71            // 消去法適用
72            for (int k = i + 1; k < n; k++)
73            {
74                double factor = matrix[k, i] / matrix[i, i];
75                for (int j = i; j < n; j++)
76                    matrix[k, j] -= factor * matrix[i, j];
77                vector[k] -= factor * vector[i];
78            }
79        }
80
81        // 後退代入で解を求める
82        double[] result = new double[n];
83        for (int i = n - 1; i >= 0; i--)
84        {
85            result[i] = vector[i];
86            for (int j = i + 1; j < n; j++)
87                result[i] -= matrix[i, j] * result[j];
88            result[i] /= matrix[i, i];
89        }
90
91        return result;
92    }
93}
94

発生している問題・分からないこと

係数aがプラスの解が欲しいのに係数aがマイナスになってしまう。

該当のソースコード

using System;

class LeastSquaresQuadratic
{
    static void Main()
    {
        // サンプルデータ（x, y）を定義
        double[] xData = { -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 };
        double[] yData = { 16, 9, 4, 1, 0, 1, 4, 9, 16, 25 };

        // 必要な各種合計値を初期化
        int n = xData.Length;
        double sumX = 0, sumY = 0, sumX2 = 0, sumX3 = 0, sumX4 = 0;
        double sumXY = 0, sumX2Y = 0;

        // 各値を計算
        for (int i = 0; i < n; i++)
        {
            double x = xData[i];
            double y = yData[i];
            sumX += x;
            sumY += y;
            sumX2 += x * x;
            sumX3 += x * x * x;
            sumX4 += x * x * x * x;
            sumXY += x * y;
            sumX2Y += x * x * y;
        }

        // 行列の各要素を計算
        double[,] matrix = {
            { sumX4, sumX3, sumX2 },
            { sumX3, sumX2, sumX },
            { sumX2, sumX, n }
        };
        double[] vector = { sumX2Y, sumXY, sumY };

        // 連立方程式を解く
        double[] coefficients = GaussianElimination(matrix, vector);

        // 結果出力
        Console.WriteLine($"y = {coefficients[0]:F3}x^2 + {coefficients[1]:F3}x + {coefficients[2]:F3}");
    }

    // ガウス消去法で連立方程式を解くメソッド
    static double[] GaussianElimination(double[,] matrix, double[] vector)
    {
        int n = vector.Length;

        for (int i = 0; i < n; i++)
        {
            // 対角要素がゼロでないようにするためのピボット選択
            int maxRow = i;
            for (int k = i + 1; k < n; k++)
            {
                if (Math.Abs(matrix[k, i]) > Math.Abs(matrix[maxRow, i]))
                    maxRow = k;
            }

            // 行の入れ替え
            for (int k = i; k < n; k++)
            {
                double tmp = matrix[maxRow, k];
                matrix[maxRow, k] = matrix[i, k];
                matrix[i, k] = tmp;
            }
            double tmp2 = vector[maxRow];
            vector[maxRow] = vector[i];
            vector[i] = tmp2;

            // 消去法適用
            for (int k = i + 1; k < n; k++)
            {
                double factor = matrix[k, i] / matrix[i, i];
                for (int j = i; j < n; j++)
                    matrix[k, j] -= factor * matrix[i, j];
                vector[k] -= factor * vector[i];
            }
        }

        // 後退代入で解を求める
        double[] result = new double[n];
        for (int i = n - 1; i >= 0; i--)
        {
            result[i] = vector[i];
            for (int j = i + 1; j < n; j++)
                result[i] -= matrix[i, j] * result[j];
            result[i] /= matrix[i, i];
        }

        return result;
    }
}

試したこと・調べたこと

teratailやGoogle等で検索した
ソースコードを自分なりに変更した
知人に聞いた
その他

上記の詳細・結果

思った通りに動かない

補足

特になし

kikukiku

2024/10/29 23:59

下記の計算方法の方が単純で簡単ではないでしょうか？こちらで実装してみてはいかがでしょうか？ https://sci-pursuit.com/math/statistics/least-square-method.html

jmdajmw

2024/10/30 01:00

やっぱ二次関数に近似する時点で無理がありますよね、、、

kikukiku

2024/10/30 01:23

２次関数だったのですね。間違えました。取り下げます。

jmdajmw

2024/10/30 01:33

いえ、二次関数に近い曲線ではありますが求めたい正解は近似値ではなくズバリぴったり合う値なので一次関数に近似してそれ以外のところでいろいろ工夫した方がよさそうな感じです。

bsdfan

2024/10/30 01:41

二次関数でも問題なく最小二乗法できるはずです。 https://sciotein.com/2022/11/01/least-squares-method-matrix/ ぱっと見た感じ作っている行列はあってそうなので、GaussianElimination に原因があるんじゃないでしょうか。GaussianElimination が正しく動いているか確かめてみてはどうでしょう。

kikukiku

2024/10/30 01:44

下記に２次関数の場合のパラメータの求め方が記述されているので、こちらで計算してみてはだめでしょうか？ http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/tech32.html

winterboum

2024/10/30 02:03

どんな値になったのかな。。ダブルだし有効桁数少ないし、きれいな放物線だし、、、桁落ちが原因とは脊椎反射では「なさそ」では有るのですが、実際系では起きる可能性高まるのでちょっとコメント。こういう場合ｘ，ｙ共に平均値出してその差分で計算するのをお薦めします。求めた解に平均値を足したのが解

jmdajmw

2024/10/30 03:11

山の頂点がズレると実データから次の目標を求めても全然狙った値にならないので最初に山の頂点を測定するひと手間が必要になります。

kikukiku

2024/10/31 06:58

元のソースの原因を追究しました。下記のソースに入れ替えたところ、誤差なく係数を求められたため、演算による誤差だと思います。 // 消去法適用 for (int k = i + 1; k < n; k++) { double factor = matrix[k, i] / matrix[i, i]; for (int j = i; j < n; j++) { matrix[k, j] -= factor * matrix[i, j]; //●計算誤差で０にならないため、極小さい数字の場合には強制的に０にする if (matrix[k, j] < 0.000001) { matrix[k, j] = 0; } } vector[k] -= factor * vector[i]; //●計算誤差で０にならないため、極小さい数字の場合には強制的に０にする if (vector[k] < 0.000001) { vector[k] = 0; } }

行動規範の内容に同意します

回答1件

ベストアンサー

下記の計算式をプログラムすると、ぴったり理論値になりました。
http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/tech32.html

C#
1        private class Point
2        {
3            public double x { get; set; }
4            public double y { get; set; }
5        }
6        private void test02()
7        {
8            var datalist = new List<Point>() {
9                new Point { x = -4, y = 16}
10                , new Point { x = -3, y = 9}
11                , new Point { x = -2, y = 4}
12                , new Point { x = -1, y = 1}
13                , new Point { x = 0, y = 0}
14                , new Point { x = 1, y = 1}
15                , new Point { x = 2, y = 4}
16                , new Point { x = 3, y = 9}
17                , new Point { x = 4, y = 16}
18                , new Point { x = 5, y = 25}
19            };
20
21            var sumX0 = datalist.Sum(p => Math.Pow(p.x, 0));
22            var sumX1 = datalist.Sum(p => Math.Pow(p.x, 1));
23            var sumX2 = datalist.Sum(p => Math.Pow(p.x, 2));
24            var sumX3 = datalist.Sum(p => Math.Pow(p.x, 3));
25            var sumX4 = datalist.Sum(p => Math.Pow(p.x, 4));
26
27            var sumY = datalist.Sum(p => p.y);
28            var sumXY = datalist.Sum(p => p.x * p.y);
29            var sumXXY = datalist.Sum(p => p.x * p.x * p.y);
30
31            var a = 
32                (
33                (sumX0 * sumX2 * sumXXY)
34                - (sumX1 * sumX1 * sumXXY)
35                + (sumX1 * sumX2 * sumXY)
36                - (sumX0 * sumX3 * sumXY)
37                + (sumX1 * sumX3 * sumY)
38                - (sumX2 * sumX2 * sumY)
39                ) / (
40                (2 * sumX1 * sumX2 * sumX3)
41                + (sumX0 * sumX2 * sumX4)
42                - (sumX1 * sumX1 * sumX4)
43                - (sumX0 * sumX3 * sumX3)
44                - (sumX2 * sumX2 * sumX2)
45                )
46                ;
47
48            var b =
49                (
50                (sumX1 * sumX2 * sumXXY)
51                - (sumX0 * sumX3 * sumXXY)
52                + (sumX0 * sumX4 * sumXY)
53                - (sumX2 * sumX2 * sumXY)
54                + (sumX2 * sumX3 * sumY)
55                - (sumX1 * sumX4 * sumY)
56                ) / (
57                (2 * sumX1 * sumX2 * sumX3)
58                + (sumX0 * sumX2 * sumX4)
59                - (sumX1 * sumX1 * sumX4)
60                - (sumX0 * sumX3 * sumX3)
61                - (sumX2 * sumX2 * sumX2)
62                )
63                ;
64
65            var c =
66                (
67                - (sumX2 * sumX2 * sumXXY)
68                + (sumX1 * sumX3 * sumXXY)
69                - (sumX1 * sumX4 * sumXY)
70                + (sumX2 * sumX3 * sumXY)
71                - (sumX3 * sumX3 * sumY)
72                + (sumX2 * sumX4 * sumY)
73                ) / (
74                (2 * sumX1 * sumX2 * sumX3)
75                + (sumX0 * sumX2 * sumX4)
76                - (sumX1 * sumX1 * sumX4)
77                - (sumX0 * sumX3 * sumX3)
78                - (sumX2 * sumX2 * sumX2)
79                )
80                ;
81
82            Console.WriteLine($"a = {a}");
83            Console.WriteLine($"b = {b}");
84            Console.WriteLine($"c = {c}");
85
86            chart1.Titles.Add("Test Chart Control");
87            chart1.ChartAreas[0].AxisX.Interval = 1;
88            chart1.ChartAreas[0].AxisY.Interval = 1;
89            chart1.Series["Series1"].ChartType = SeriesChartType.Point;
90            chart1.Series["Series1"].Color = Color.DarkBlue;
91
92            chart1.Series.Add("2次曲線");
93            chart1.Series["2次曲線"].ChartType = SeriesChartType.Line;
94            chart1.Series["2次曲線"].Color = Color.Aqua;
95
96            foreach (var p in datalist)
97            {
98                chart1.Series["Series1"].Points.AddXY(p.x, p.y);
99                chart1.Series["2次曲線"].Points.AddXY(p.x, a * p.x * p.x + b * p.x + c);
100            }
101        }