Keyword: 最小二乗, 双３次, スプライン, フィット
概要
本サンプルは最小二乗双３次スプライン曲面フィットを行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて最小二乗双３次スプライン曲面フィットを行い、スプラインの値を求めて出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_2d_spline_fit_panel() のExampleコードです。本サンプル及び関数の詳細情報は nag_2d_spline_fit_panel のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_2d_spline_fit_panel のマニュアルページを参照）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
| このデータをダウンロード |
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 のマニュアルページを参照）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
| この出力例をダウンロード |
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等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
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 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250
| このソースコードをダウンロード |
/* 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;
}