Keyword: 最小二乗, 双３次, スプライン, フィット
概要
本サンプルは最小二乗双３次スプライン曲面フィットを行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて最小二乗双３次スプライン曲面フィットを行い、スプラインの値を求めて出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_2d_spline_fit_panel() のExampleコードです。本サンプル及び関数の詳細情報は nag_2d_spline_fit_panel のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_2d_spline_fit_panel のマニュアルページを参照）| このデータをダウンロード |
nag_2d_spline_fit_panel (e02dac) Example Program Data
0.000001
30
8
10
-0.52 0.60 0.93 10.
-0.61 -0.95 -1.79 10.
0.93 0.87 0.36 10.
0.09 0.84 0.52 10.
0.88 0.17 0.49 10.
-0.70 -0.87 -1.76 10.
1.00 1.00 0.33 1.
1.00 0.10 0.48 1.
0.30 0.24 0.65 1.
-0.77 -0.77 -1.82 1.
-0.23 0.32 0.92 1.
-1.00 1.00 1.00 1.
-0.26 -0.63 8.88 1.
-0.83 -0.66 -2.01 1.
0.22 0.93 0.47 1.
0.89 0.15 0.49 1.
-0.80 0.99 0.84 1.
-0.88 -0.54 -2.42 1.
0.68 0.44 0.47 1.
-0.14 -0.72 7.15 1.
0.67 0.63 0.44 1.
-0.90 -0.40 -3.34 1.
-0.84 0.20 2.78 1.
0.84 0.43 0.44 1.
0.15 0.28 0.70 1.
-0.91 -0.24 -6.52 1.
-0.35 0.86 0.66 1.
-0.16 -0.41 2.32 1.
-0.35 -0.05 1.66 1.
-1.00 -1.00 -1.00 1.
-0.5
0.0
- １行目はタイトル行で読み飛ばされます。
- ２行目に一次方程式の有効ランクを決定するための閾値(eps)を指定しています。
- ３行目にデータ点の数(m)を指定しています。
- ４行目に変数xのノット数の合計(px)を指定しています。
- ５行目に変数yのノット数の合計(py)を指定しています。
- ６〜３５行目にデータ点yの座標、xの座標、fの座標、重み(w)を指定しています。
- ３６〜３７行目に変数yに関する内部ノット(mu)を指定しています。
出力結果
（本関数の詳細はnag_2d_spline_fit_panel のマニュアルページを参照）| この出力例をダウンロード |
nag_2d_spline_fit_panel (e02dac) Example Program Results
Interior y -knots
-0.5000
0.0000
Interior x -knots
None
Sum of squares of residual RHS 1.467e+01
Rank 22
x and y have been interchanged
X Y Data Fit Residual
-0.9500 -0.6100 -1.7900 -1.7931 -3.126e-03
-0.8700 -0.7000 -1.7600 -1.7521 7.893e-03
-0.7700 -0.7700 -1.8200 -2.4301 -6.101e-01
-0.6300 -0.2600 8.8800 7.6346 -1.245e+00
-0.6600 -0.8300 -2.0100 -1.5815 4.285e-01
-0.5400 -0.8800 -2.4200 -2.6795 -2.595e-01
-0.7200 -0.1400 7.1500 7.5708 4.208e-01
-1.0000 -1.0000 -1.0000 -1.0228 -2.277e-02
-0.4000 -0.9000 -3.3400 -4.6955 -1.356e+00
-0.2400 -0.9100 -6.5200 -4.7072 1.813e+00
-0.4100 -0.1600 2.3200 2.7039 3.839e-01
-0.0500 -0.3500 1.6600 2.2865 6.265e-01
0.6000 -0.5200 0.9300 0.9441 1.407e-02
0.8700 0.9300 0.3600 0.3529 -7.097e-03
0.8400 0.0900 0.5200 0.5024 -1.761e-02
0.1700 0.8800 0.4900 0.4705 -1.951e-02
1.0000 1.0000 0.3300 0.6315 3.015e-01
0.1000 1.0000 0.4800 1.4910 1.011e+00
0.2400 0.3000 0.6500 0.9241 2.741e-01
0.3200 -0.2300 0.9200 -0.3692 -1.289e+00
1.0000 -1.0000 1.0000 1.0835 8.352e-02
0.9300 0.2200 0.4700 1.4912 1.021e+00
0.1500 0.8900 0.4900 0.4414 -4.861e-02
0.9900 -0.8000 0.8400 0.5495 -2.905e-01
0.4400 0.6800 0.4700 1.5862 1.116e+00
0.6300 0.6700 0.4400 0.6288 1.888e-01
0.2000 -0.8400 2.7800 1.7123 -1.068e+00
0.4300 0.8400 0.4400 0.6888 2.488e-01
0.2800 0.1500 0.7000 0.7713 7.134e-02
0.8600 -0.3500 0.6600 0.9347 2.747e-01
Sum of squared residuals 1.467e+01
Spline coefficients
-1.0228 115.4668 -433.5558 -68.1973
24.8426 -140.1485 258.5042 15.6756
-29.4878 132.2933 -173.5103 20.0983
9.9575 -51.6200 67.6666 -5.8765
10.0577 4.7543 -15.3533 -0.3260
1.0835 -2.7932 7.7708 0.6315
- ３〜５行目に変数yの内部ノットが出力されています。
- ７〜８行目に変数xの内部ノットがないことが示されています。
- １０行目に残差平方和が出力されています。
- １２行目にランクが出力されています。
- １４行目にxとyが入れ替えられたことが示されています。
- １６〜４６行目にx、y、f、フィットされた値、残差が出力されています。
- ４８行目に残差平方和が出力されています。
- ５０〜５６行目にスプライン係数が出力されています。
ソースコード
（本関数の詳細はnag_2d_spline_fit_panel のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_2d_spline_fit_panel (e02dac) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage02.h>
int main(void)
{
/* Initialized data */
char label[] = "xy";
/* Scalars */
double d, eps, sigma, sum, temp;
Integer exit_status = 0, i, iadres, itemp, j, m, nc, np, npoint, px, py,
rank;
/* Arrays */
double *dl = 0, *f = 0, *ff = 0, *lamda = 0, *mu = 0, *w = 0, *x = 0;
double *y = 0;
Integer *point = 0;
/* Nag Types */
Nag_2dSpline spline;
NagError fail;
exit_status = 0;
INIT_FAIL(fail);
/* Initialize spline */
spline.lamda = 0;
spline.mu = 0;
spline.c = 0;
printf("nag_2d_spline_fit_panel (e02dac) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
while (scanf("%lf", &eps) != EOF && exit_status == 0)
{
/* Read data, interchanging X and Y axes if PX.LT.PY */
scanf("%ld%*[^\n] ", &m);
if (m > 1) {
/* Allocate memory */
if (!(f = nAG_ALLOC(m, double)) ||
!(ff = nAG_ALLOC(m, double)) ||
!(w = 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%*[^\n] ", &px, &py);
if (px < 8 && py < 8) {
printf("px or py is too small.\n");
exit_status = 1;
goto END;
}
nc = (px - 4) * (py - 4);
np = (px - 7) * (py - 7);
npoint = m + (px - 7) * (py - 7);
/* Allocate memory */
if (!(dl = nAG_ALLOC(nc, double)) ||
!(point = nAG_ALLOC(npoint, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
if (px < py) {
itemp = px;
px = py;
py = itemp;
itemp = 1;
/* Allocate memory */
if (!(lamda = nAG_ALLOC(px, double)) || !(mu = nAG_ALLOC(py, double))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 0; i < m; ++i)
scanf("%lf%lf%lf%lf", &y[i], &x[i], &f[i], &w[i]);
scanf("%*[^\n] ");
if (py > 8) {
for (j = 4; j < py - 4; ++j)
scanf("%lf", &mu[j]);
scanf("%*[^\n] ");
}
if (px > 8) {
for (j = 4; j < px - 4; ++j)
scanf("%lf", &lamda[j]);
scanf("%*[^\n] ");
}
}
else {
/* Allocate memory */
if (!(lamda = nAG_ALLOC(px, double)) || !(mu = nAG_ALLOC(py, double))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
itemp = 0;
for (i = 0; i < m; ++i)
scanf("%lf%lf%lf%lf", &x[i], &y[i], &f[i], &w[i]);
scanf("%*[^\n] ");
if (px > 8) {
for (j = 4; j < px - 4; ++j)
scanf("%lf", &lamda[j]);
scanf("%*[^\n] ");
}
if (py > 8) {
for (j = 4; j < py - 4; ++j)
scanf("%lf", &mu[j]);
scanf("%*[^\n] ");
}
}
printf("\nInterior %1.1s -knots\n", label + itemp);
for (j = 4; j < px - 4; ++j)
printf("%11.4f\n", lamda[j]);
if (px == 8)
printf("None\n");
printf("\nInterior %1.1s -knots\n", label + (2 - itemp - 1));
for (j = 4; j < py - 4; ++j)
printf("%1s%11.4f\n", "", mu[j]);
if (py == 8)
printf("None\n");
/* nag_2d_panel_sort (e02zac).
* Sort two-dimensional data into panels for fitting bicubic
* splines
*/
nag_2d_panel_sort(px, py, lamda, mu, m, x, y, point, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_panel_sort (e02zac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Fit bicubic spline to data points */
spline.nx = px;
spline.ny = py;
if (!(spline.c = nAG_ALLOC((spline.nx - 4) * (spline.ny - 4), double)) ||
!(spline.lamda = nAG_ALLOC(spline.nx, double)) ||
!(spline.mu = nAG_ALLOC(spline.ny, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 0; i < spline.nx; i++)
spline.lamda[i] = lamda[i];
for (i = 0; i < spline.ny; i++)
spline.mu[i] = mu[i];
/* nag_2d_spline_fit_panel (e02dac).
* Least squares surface fit, bicubic splines
*/
nag_2d_spline_fit_panel(m, x, y, f, w, point, dl, eps, &sigma, &rank,
&spline, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_fit_panel (e02dac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("\nSum of squares of residual RHS%16.3e\n", sigma);
printf("\nRank%5ld\n", rank);
/* nag_2d_spline_eval (e02dec).
* Evaluation of bicubic spline, at a set of points
*/
nag_2d_spline_eval(m, x, y, ff, &spline, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_eval (e02dec).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
sum = 0.;
if (itemp == 1)
printf("\nx and y have been interchanged\n\n");
/*Output data points, fitted values and residuals */
printf(" X Y Data Fit Residual\n");
for (i = 0; i < np; ++i) {
iadres = i + m;
while ((iadres = point[iadres] - 1) >= 0) {
temp = ff[iadres] - f[iadres];
printf("%11.4f%11.4f%11.4f%11.4f%12.3e\n", x[iadres],
y[iadres], f[iadres], ff[iadres], temp);
/* Computing 2nd power */
d = temp * w[iadres];
sum += d * d;
}
}
printf("\nSum of squared residuals%16.3e\n", sum);
printf("\nSpline coefficients\n");
for (i = 0; i < px - 4; ++i) {
for (j = 0; j < py - 4; ++j)
printf("%11.4f", spline.c[i * (py - 4) + j]);
printf("\n");
}
END:
nAG_FREE(dl);
nAG_FREE(f);
nAG_FREE(ff);
nAG_FREE(lamda);
nAG_FREE(mu);
nAG_FREE(w);
nAG_FREE(x);
nAG_FREE(y);
nAG_FREE(point);
nAG_FREE(spline.lamda);
nAG_FREE(spline.mu);
nAG_FREE(spline.c);
}
return exit_status;
}