Karinto/CurveFitter.cs

   1 /*\r
   2  *      Karinto Library Project\r
   3  *\r
   4  *      This software is distributed under a zlib-style license.\r
   5  *      See license.txt for more information.\r
   6  */\r
   7 \r
   8 using System;\r
   9 using System.Collections.Generic;\r
  10 using System.Diagnostics;\r
  11 \r
  12 namespace Karinto\r
  13 {\r
  14     /// <summary>\r
  15     ///     最小二乗法による近似多項式を求める\r
  16     /// </summary>\r
  17     public static class CurveFitter\r
  18     {\r
  19         public static Polynomial Fit(PointList target, int order)\r
  20         {\r
  21             switch (order)\r
  22             {\r
  23                 case 2:\r
  24                     return FitQuadratic(target);\r
  25                 case 1:\r
  26                     return FitLinear(target);\r
  27             }\r
  28             throw new NotImplementedException();\r
  29         }\r
  30 \r
  31         public static Polynomial FitLinear(PointList target)\r
  32         {\r
  33             // an = Σ x[i]^n\r
  34             // bn = Σ x[i]^n * y[i]\r
  35             double a2 = 0, a1 = 0, a0 = target.Count;\r
  36             double b1 = 0, b0 = 0;\r
  37             foreach (Point p in target)\r
  38             {\r
  39                 double xi = p.X;\r
  40                 double yi = p.Y;\r
  41                 b0 += yi;\r
  42                 b1 += yi * xi;\r
  43                 a1 += xi;\r
  44                 a2 += xi * xi;\r
  45             }\r
  46             // ( a2 a1 )・( k1 ) = ( b1 )\r
  47             // ( a1 a0 )  ( k0 )   ( b0 )\r
  48             double det = a0 * a2 - a1 * a1; // 係数行列の行列式\r
  49             if (Math.Abs(det) < Double.Epsilon)\r
  50             {\r
  51                 // 不定の場合はどうしよう\r
  52                 throw new NotImplementedException();\r
  53             }\r
  54             double k1 = a0 * b1 - a1 * b0;\r
  55             double k0 = -a1 * b1 + a2 * b0;\r
  56             double[] c = new Double[] { k0 / det, k1 / det };\r
  57             return new Polynomial(c);\r
  58         }\r
  59 \r
  60         public static Polynomial FitQuadratic(PointList target)\r
  61         {\r
  62             // an = Σ x[i]^n\r
  63             // bn = Σ x[i]^n * y[i]\r
  64             double a4 = 0, a3 = 0, a2 = 0, a1 = 0, a0 = target.Count;\r
  65             double b2 = 0, b1 = 0, b0 = 0;\r
  66             foreach (Point p in target)\r
  67             {\r
  68                 double xi = p.X;\r
  69                 double yi = p.Y;\r
  70                 double xx = xi * xi;\r
  71                 double yx = yi * xi;\r
  72                 b0 += yi;\r
  73                 b1 += yx;\r
  74                 b2 += yx * xi;\r
  75                 a1 += xi;\r
  76                 a2 += xx;\r
  77                 a3 += xx * xi;\r
  78                 a4 += xx * xx;\r
  79             }\r
  80             // 3元連立一次方程式の解をクラメルの公式で考える\r
  81             // ( a4 a3 a2 )  ( k2 )   ( b2 )\r
  82             // | a3 a2 a1 |・| k1 | = | b1 |\r
  83             // ( a2 a1 a0 )  ( k0 )   ( b0 )\r
  84             double det = (a0 * a2 - a1 * a1) * a4 - a0 * a3 * a3 + 2 * a1 * a2 * a3 - a2 * a2 * a2; // 係数行列の行列式\r
  85             if (Math.Abs(det) < Double.Epsilon)\r
  86             {\r
  87                 // 不定の場合は直線フィッティング\r
  88                 return FitLinear(target);\r
  89             }\r
  90             double k2 = (a0 * a2 - a1 * a1) * b2 + (a1 * a2 - a0 * a3) * b1 + (a1 * a3 - a2 * a2) * b0;\r
  91             double k1 = (a1 * a2 - a0 * a3) * b2 + (a0 * a4 - a2 * a2) * b1 + (a2 * a3 - a1 * a4) * b0;\r
  92             double k0 = (a1 * a3 - a2 * a2) * b2 + (a2 * a3 - a1 * a4) * b1 + (a2 * a4 - a3 * a3) * b0;\r
  93             double[] c = new double[] { k0 / det, k1 / det, k2 / det };\r
  94             return new Polynomial(c);\r
  95         }\r
  96 \r
  97         public static double CalucurateR2(PointList target, Polynomial approx)\r
  98         {\r
  99             double[] yArray = approx.GetValues(target.XArray);\r
 100             PointList yy = new PointList(target.YArray, yArray);\r
 101             return yy.R2;\r
 102         }\r
 103     }\r
 104 }\r