最小二乗3次スプライン曲線フィット

C言語によるサンプルソースコード : 使用関数名:nag_1d_spline_fit_knots (e02bac)

ホーム > 製品 > nAG数値計算ライブラリ > サンプルソースコード集 > 最小二乗3次スプライン曲線フィット (C言語/C++)

Keyword: 最小二乗, 3次, スプライン, フィット

概要

本サンプルは最小二乗3次スプライン曲線フィットを行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて最小二乗3次スプライン曲線フィットを行い、スプラインの値を求めて出力します。

3次スプラインのデータ 

※本サンプルはnAG Cライブラリに含まれる関数 nag_1d_spline_fit_knots() のExampleコードです。本サンプル及び関数の詳細情報は nag_1d_spline_fit_knots のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

(本関数の詳細はnag_1d_spline_fit_knots のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

このデータをダウンロード
nag_1d_spline_fit_knots (e02bac) Example Program Data
  14
   5   2
     1.50
     2.60
     4.00
     8.00
     0.20     0.00     0.20
     0.47     2.00     0.20
     0.74     4.00     0.30
     1.09     6.00     0.70
     1.60     8.00     0.90
     1.90     8.62     1.00
     2.60     9.10     1.00
     3.10     8.90     1.00
     4.00     8.15     0.80
     5.15     7.00     0.50
     6.17     6.00     0.70
     8.00     4.54     1.00
    10.00     3.39     1.00
    12.00     2.56     1.00

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目にデータ点の数(m)を指定しています。
  • 3行目にスプラインの区間の数(ncap)、重みづけをするかどうかのフラグ(wght)を指定しています。"1"は重みづけをしないことを意味し、"1"以外は重みづけをすることを意味します。
  • 4~7行目にノット(lamda)を指定しています。
  • 8~21行目に独立変数xの値、従属変数yの値、重み(weights)を指定しています。

出力結果

(本関数の詳細はnag_1d_spline_fit_knots のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

この出力例をダウンロード
nag_1d_spline_fit_knots (e02bac) Example Program Results

Number of distinct knots = 6

Distinct knots located at 

  0.2000   1.5000   2.6000   4.0000   8.0000  12.0000


    J    B-spline coeff c

    1        -0.0465
    2         3.6150
    3         8.5724
    4         9.4261
    5         7.2716
    6         4.1207
    7         3.0822
    8         2.5597

Residual sum of squares =    1.78e-03

Cubic spline approximation and residuals
  r        Abscissa        Weight       Ordinate          Spline      Residual

  1         0.2000         0.2000         0.0000         -0.0465      -4.7e-02
            0.3350                                        1.0622
  2         0.4700         0.2000         2.0000          2.1057       1.1e-01
            0.6050                                        3.0817
  3         0.7400         0.3000         4.0000          3.9880      -1.2e-02
            0.9150                                        5.0558
  4         1.0900         0.7000         6.0000          5.9983      -1.7e-03
            1.3450                                        7.1376
  5         1.6000         0.9000         8.0000          7.9872      -1.3e-02
            1.7500                                        8.3544
  6         1.9000         1.0000         8.6200          8.6348       1.5e-02
            2.2500                                        9.0076
  7         2.6000         1.0000         9.1000          9.0896      -1.0e-02
            2.8500                                        9.0353
  8         3.1000         1.0000         8.9000          8.9125       1.2e-02
            3.5500                                        8.5660
  9         4.0000         0.8000         8.1500          8.1321      -1.8e-02
            4.5750                                        7.5592
 10         5.1500         0.5000         7.0000          6.9925      -7.5e-03
            5.6600                                        6.5010
 11         6.1700         0.7000         6.0000          6.0255       2.6e-02
            7.0850                                        5.2292
 12         8.0000         1.0000         4.5400          4.5315      -8.5e-03
            9.0000                                        3.9045
 13        10.0000         1.0000         3.3900          3.3928       2.8e-03
           11.0000                                        2.9574
 14        12.0000         1.0000         2.5600          2.5597      -3.5e-04

  • 3行目にノットの数が出力されています。
  • 5~7行目にノットの位置が出力されています。
  • 10~19行目にBスプライン係数が出力されています。
  • 21行目に残差平方和が出力されています。
  • 23~52行目にインデックス、独立変数xの値、重み、従属変数yの値、3次スプラインフィットの値、、残差が出力されています。

ソースコード

(本関数の詳細はnag_1d_spline_fit_knots のマニュアルページを参照)

※本サンプルソースコードはnAG数値計算ライブラリ(Windows, Linux, MAC等に対応)の関数を呼び出します。
サンプルのコンパイル及び実行方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

このソースコードをダウンロード
/* nag_1d_spline_fit_knots (e02bac) Example Program.
 *
 * CLL6I261D/CLL6I261DL Version.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nage02.h>

int main(void)
{
  Integer exit_status = 0, j, m, ncap, ncap7, r, wght;
  NagError fail;
  Nag_Spline spline;
  double fit, ss, *weights = 0, *x = 0, xarg, *y = 0;

  INIT_FAIL(fail);

  /* Initialize spline */
  spline.lamda = 0;
  spline.c = 0;

  printf("nag_1d_spline_fit_knots (e02bac) Example Program Results\n");
  scanf("%*[^\n]"); /* Skip heading in data file */
  while (scanf("%ld", &m) != EOF)
  {
    if (m >= 4) {
      if (!(weights = nAG_ALLOC(m, double)) ||
          !(x = nAG_ALLOC(m, double)) || !(y = nAG_ALLOC(m, double))
             )
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    }
    else {
      printf("Invalid m.\n");
      exit_status = 1;
      goto END;
    }
    scanf("%ld%ld", &ncap, &wght);
    if (ncap > 0) {
      ncap7 = ncap + 7;
      spline.n = ncap7;
      if (!(spline.lamda = nAG_ALLOC(ncap7, double)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    }
    else {
      printf("Invalid ncap.\n");
      exit_status = 1;
      goto END;
    }
    for (j = 4; j < ncap + 3; ++j)
      scanf("%lf", &(spline.lamda[j]));
    for (r = 0; r < m; ++r) {
      if (wght == 1) {
        scanf("%lf%lf", &x[r], &y[r]);
        weights[r] = 1.0;
      }
      else
        scanf("%lf%lf%lf", &x[r], &y[r], &weights[r]);
    }
    /* nag_1d_spline_fit_knots (e02bac).
     * Least squares curve cubic spline fit (including
     * interpolation), one variable
     */
    nag_1d_spline_fit_knots(m, x, y, weights, &ss, &spline, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_1d_spline_fit_knots (e02bac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    printf("\nNumber of distinct knots = %ld\n\n", ncap + 1);
    printf("Distinct knots located at \n\n");
    for (j = 3; j < ncap + 4; j++)
      printf("%8.4f%s", spline.lamda[j],
             (j - 3) % 6 == 5 || j == ncap + 3 ? "\n" : " ");
    printf("\n\n    J    B-spline coeff c\n\n");
    for (j = 0; j < ncap + 3; ++j)
      printf("    %ld  %13.4f\n", j + 1, spline.c[j]);
    printf("\nResidual sum of squares = ");
    printf("%11.2e\n\n", ss);
    printf("Cubic spline approximation and residuals\n");
    printf("  r        Abscissa        Weight       Ordinate"
           "          Spline      Residual\n\n");
    for (r = 0; r < m; ++r) {
      /* nag_1d_spline_evaluate (e02bbc).
       * Evaluation of fitted cubic spline, function only
       */
      nag_1d_spline_evaluate(x[r], &fit, &spline, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_1d_spline_evaluate (e02bbc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }

      printf("%3ld    %11.4f    %11.4f    %11.4f     %11.4f"
             "    %10.1e\n", r + 1, x[r], weights[r], y[r], fit, fit - y[r]);
      if (r < m - 1) {
        xarg = (x[r] + x[r + 1]) * 0.5;
        /* nag_1d_spline_evaluate (e02bbc), see above. */
        nag_1d_spline_evaluate(xarg, &fit, &spline, &fail);
        if (fail.code != NE_NOERROR) {
          printf("Error from nag_1d_spline_evaluate (e02bbc).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }
        printf("    %14.4f             %33.4f\n", xarg, fit);
      }
    }
    /* Free memory used by spline */
    nAG_FREE(spline.lamda);
    nAG_FREE(spline.c);
  END:
    nAG_FREE(weights);
    nAG_FREE(x);
    nAG_FREE(y);
  }
  return exit_status;
}


関連情報
© 日本ニューメリカルアルゴリズムズグループ株式会社 2024
Privacy Policy  /  Trademarks