fork(1) download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. void gauss3(double A[3][3], double b[3], double x[3]) {
  5.     for (int i = 0; i < 3; i++) {
  6.         // ピボット行の正規化
  7.         double pivot = A[i][i];
  8.         for (int j = 0; j < 3; j++) A[i][j] /= pivot;
  9.         b[i] /= pivot;
  10.  
  11.         // 他の行の消去 (掃き出し処理)
  12.         for (int k = 0; k < 3; k++) {
  13.             if (k == i) continue;
  14.             double factor = A[k][i];
  15.             for (int j = 0; j < 3; j++) A[k][j] -= factor * A[i][j];
  16.             b[k] -= factor * b[i];
  17.         }
  18.     }
  19.  
  20.     for (int i = 0; i < 3; i++) x[i] = b[i];
  21. }
  22.  
  23. int main() {
  24.     // 与えられたデータセット (8組)
  25.     double xval[8] = {0.15708, 0.23982, 0.37400, 0.57120, 0.82674, 1.04720, 1.23200, 1.43452};
  26.     double yval[8] = {0.98769, 0.97138, 0.93087, 0.84125, 0.67728, 0.50000, 0.33236, 0.13586};
  27.  
  28.     double S2 = 0, S4 = 0, S6 = 0, S8 = 0;
  29.     double Sy = 0, S2y = 0, S4y = 0;
  30.  
  31.     // 正規方程式に必要な各合計値の計算
  32.     for (int i = 0; i < 8; i++) {
  33.         double x2 = xval[i] * xval[i];
  34.         double x4 = x2 * x2;
  35.         double x6 = x4 * x2;
  36.         double x8 = x4 * x4;
  37.  
  38.         S2 += x2;
  39.         S4 += x4;
  40.         S6 += x6;
  41.         S8 += x8;
  42.  
  43.         Sy += yval[i];
  44.         S2y += x2 * yval[i];
  45.         S4y += x4 * yval[i];
  46.     }
  47.  
  48.     // 係数行列 A と右辺ベクトル b の設定
  49.     double A[3][3] = {
  50.         {8,  S2, S4},
  51.         {S2, S4, S6},
  52.         {S4, S6, S8}
  53.     };
  54.  
  55.     double b[3] = {Sy, S2y, S4y};
  56.     double sol[3];
  57.  
  58.     // 掃き出し法 (ガウスの消去法) で連立方程式を解く
  59.     gauss3(A, b, sol);
  60.  
  61.     // 結果の出力
  62.     printf("=========================================\n");
  63.     printf("モデル (2): y = b1 + b2*x^2 + b3*x^4 計算結果\n");
  64.     printf("=========================================\n");
  65.     printf("求められた係数:\n");
  66.     printf("  b1 = %.6f\n", sol[0]);
  67.     printf("  b2 = %.6f\n", sol[1]);
  68.     printf("  b3 = %.6f\n", sol[2]);
  69.  
  70.     // 残差平方和 (誤差の2乗和 E) の計算
  71.     double E = 0;
  72.     for (int i = 0; i < 8; i++) {
  73.         double x2 = xval[i] * xval[i];
  74.         double x4 = x2 * x2;
  75.         double yhat = sol[0] + sol[1] * x2 + sol[2] * x4;
  76.         double diff = yval[i] - yhat;
  77.         E += diff * diff;
  78.     }
  79.  
  80.     printf("\n近似式による偏差の2乗和 (誤差 E):\n");
  81.     printf("  E = %.10f\n", E);
  82.  
  83.     // Excelプロット用の予測値出力
  84.     printf("\n近似式による y の予測値 (Excelグラフ作成用):\n");
  85.     printf("x 値\t\t計算された y 予測値\n");
  86.     for (int i = 0; i < 8; i++) {
  87.         double x2 = xval[i] * xval[i];
  88.         double x4 = x2 * x2;
  89.         double yhat = sol[0] + sol[1] * x2 + sol[2] * x4;
  90.         printf("%.5f\t%.6f\n", xval[i], yhat);
  91.     }
  92.     printf("=========================================\n");
  93.  
  94.     return 0;
  95. }
Success #stdin #stdout 0s 5304KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=========================================
モデル (2): y = b1 + b2*x^2 + b3*x^4 計算結果
=========================================
求められた係数:
  b1 = 0.999699
  b2 = -0.496930
  b3 = 0.037543

近似式による偏差の2乗和 (誤差 E):
  E = 0.0000005016

近似式による y の予測値 (Excelグラフ作成用):
x 値		計算された y 予測値
0.15708	0.987460
0.23982	0.971243
0.37400	0.930925
0.57120	0.841562
0.82674	0.677587
1.04720	0.499900
1.23200	0.331937
1.43452	0.136076
=========================================
https://ideone.com/88d8sH
language:
C (gcc 8.3)
created:
16 hours ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!