|
Maple 16 の補間プロットと平滑化
このワークシートでは、2-D データと 3-D データの補間と平滑化について、次元、タスク、およびデータの規則性で分類して説明します。
共通するタスクは、データポイントを通過または近似する曲線または面として、データのプロットだけを生成することです。
もう 1 つのタスクは、データを近似し、元のデータポイントの近くにある個々のポイントの値をクエリすることができるプロシージャや区分関数スプライン式を生成することです。
もう 1 つの選択肢は、指定したデータポイントを通過する曲線または面を生成する (つまりデータポイントを補間する) か、指定したデータポイントを (平滑化によって) 単に近似するかということです。データにノイズや測定誤差が含まれているかどうかを確認してください。依存データに誤差やノイズが含まれることが予想される場合、平滑化した面や曲線をデータに適合させることが妥当です。この場合、平滑化したカーブフィットが、必ずしも、指定したデータポイントの従属データを通過したり、そのデータと一致するとは限りません。一方、データが完全に正確なことが予想される場合、指定したすべてのデータポイントを通過する、補間処理した曲線または面を生成することが妥当です。
データから直接、曲線または面をプロットすることも、または最初にプロシージャや (2 次元の場合は) 区分関数スプラインを作成してから、plot または plot3d コマンドにそのプロシージャやスプラインを提供することで、間接的にプロットを実行することもできます。作成したプロシージャは、個々の 1 次元または 2 次元独立データポイントを受け付け、その入力に対する従属スカラー値を計算します。
3D データの関連するもう 1 つの区分は、データの独立部が、2D 平面の規則的な格子上に配置されるか、双方向で不規則に配置されているポイントによって構成されているかという点です。
このワークシートでは、2D プロットと 3D プロットに対するこのような手法の詳細な説明と比較を行います。
|
>
|
restart;
with(LinearAlgebra): with(CurveFitting): with(plots): with(Statistics):
randomize():
|
|
|
1 次元の場合
|
|
X データポイントは、1 ~ 15 の範囲の整数として均一に取得されます。しかし、これらのデータポイントは順番が正しい限り、不均一な間隔に配置されるほうが望ましいこともあります。
この例では、data2D は 15 x 2 の行列で、最初の列は順序のある x 方向のデータポイント、2 番目の列はその y の値です。
|
>
|
data2D:=< <($1..15)> | RandomVector(15,generator=0.0..5.0) >:
|
線形補間を生成する plot コマンドを次に示します。このコマンドの目的は、上記コマンドの代わりにスムーズな曲線を生成することです。
|
>
|
display( plot(data2D, color=black),
pointplot(data2D, symbolsize=10, color=red) );
|
ここでは、補間スプラインの区分関数式を生成し、この式をプロットします。
|
>
|
f := Spline(data2D, v):
display(plot(f, v=1..15, color=black),
pointplot(data2D, symbolsize=10, color=red) );
|
区分関数スプラインの記号式を生成する上記手法は、最もシンプルな手法です。この手法では、プロッターは任意の数の補間ポイント数を選択できます。複数のオプションを指定して、Spline コマンドを呼び出すことができます。しかし多数の元のデータポイントを処理するには、この手法は効率的ではありません。この手法で生成される区分関数式は場合によって非常に処理しづらいものとなり、多くのポイントでの数値評価コストは大きく増加します。次に、ArrayInterpolation コマンドを使った、補間対象となるポイントの事前作成済みメッシュの使用方法について説明します。
ここでは固定メッシュの生成、対象ポイントでの補間、および生成された値のプロットを実行します。元のデータポイントが 15 個存在するため、1/10 の間隔で均一に配置されたポイントを使って、150 ポイントの 精細メッシュ A1 を作成します。
これは時間がかからず、メモリ消費量も少ない手法です。しかし、評価ポイントはX 軸に沿って均等に配置されるため、見た目が美しいプロットが生成されるとは限りません。曲線の中には、この手法で得られる結果と、デフォルトの適用 2D プロッタで得られる結果の違いが著しく大きなものもあります。
|
>
|
A1:=<(i/10$i=1..150)>:
f:=ArrayInterpolation(data2D, A1, method=spline):
g:=ArrayInterpolation(data2D, A1, method=cubic):
drawfun := f_ ->
display(plot(<A1|f_>, color=black),
pointplot(data2D, symbolsize=10, color=red)):
|
|
>
|
display(Array([drawfun(f), drawfun(g)]), view=[0..15, 0..5]);
|
|
|
補間プロシージャの取得
|
|
上記の方式では、事前作成済みベクトル A1 を使用します。これは固定された均一の補間ポイントで構成されています。また特定の数値入力に対して単一のスカラー数値出力を生成するプロシージャを指定して plot コマンドを実行した場合、補間ポイントを動的に選択することもできます。次の例では、前述のように動作する演算子 F をデモンストレーションします。
F は元のデータを使用しますが、補間ポイントを前もって提供する必要はないという点に注意してください。この手法をより効率的に実装する方法については、「付録」を参照してください。
2-D プロットの場合、Maple はデフォルトで適応プロットを実行します。この場合、変化する関数の検出方法に応じて、特定の領域でプロットする関数の評価頻度が増加します。plot コマンドの adaptive オプションの詳細については、plot,options のヘルプページを参照してください。
|
>
|
F:=x->CurveFitting:-ArrayInterpolation(data2D,
Array(1..1, 1..1, [[x]]), method=spline)[1]:
display( plot(F, 0..15, color=black),
pointplot(data2D, symbolsize=10, color=red), view=[0..15, 0..5] );
|
|
|
一般的な使用方法
|
|
その他の一般的な計算で、F などの補間プロシージャをスカラー値関数として使用することもできます。
|
>
|
Optimization:-Minimize(t->F(t), 0..15, method=branchandbound);
|
|
![[-0.459382095810952e-1, Vector(1, {(1) = 8.30778366525865})]](/support/helpjp/helpview.aspx?si=2547/file01547/math212.png)
| (1.1.1.1) |
|
>
|
fsolve(t->F(t)-2, 0..15);
|
|
| (1.1.1.2) |
|
|
|
|
|
2 次元の場合
|
|
以下に、データ補間 (指定したデータポイントと正確に一致) とデータ平滑化 (ノイズデータを近似) の違いを示します。
|
|
平滑化
|
|
このセクションでは、1 次元の依存データに誤差成分が存在するか、ノイズが含まれていることが判明している、2 次元独立データの集合があるものとします。ノイズが存在するということは、そのデータが示す面がどのようなものであっても、その面が必ずしも指定した依存データポイントを通過するとは限らないということです。
この場合、lowess アルゴリズムを使ってデータを数値的に平滑化することで、面の近似が行われます。指定したポイントをプロットする場合、データポイントの周囲に存在するウィンドウを使って、低次に重み付けられた、局所カーブフィットを計算します。
|
|
ScatterPlot3D
|
|
統計パッケージの ScatterPlot3D コマンドを実行すると、2D ノイズデータをスムーズにプロットできます。独立データは、規則的に配置する必要がなくなり、n x 3 の配列または行列として提供されます。各 n 行は、個々のポイントを表しています。列は x、y、z の各値です。
次に例を示します。最初に、x-y 平面へのデータポイントの投影について説明するため、2D 従属データ値のレイアウトと間隔を表示します。
|
>
|
X := Sample(Uniform(-50,50),175):
|
|
>
|
Y := Sample(Uniform(-50,50),175):
|
|
>
|
Zerror := Sample(Normal(0,100),175):
|
|
>
|
Z := Array(1..175,(i)->-(sin(Y[i]/20)*(X[i]-6)^2+(Y[i]-7)^2+Zerror[i])):
|
|
>
|
XYZ := Matrix([[X],[Y],[Z]],datatype=float[8])^%T;
|
|
![XYZ := Vector(4, {(1) = ` 175 x 3 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})](/support/helpjp/helpview.aspx?si=2547/file01547/math271.png)
| (2.1.1.1) |
|
>
|
ScatterPlot3D(XYZ, axes=box, orientation=[20,0,0]);
|
|
>
|
ScatterPlot3D(XYZ, lowess, grid=[25,25], axes=box, orientation=[20,70,0]);
|
次に別の例を示します。この例では、ファイルから n x 3 のデータが読み込まれます。
|
>
|
M:=ImportMatrix(cat(kernelopts(mapledir),"/data/plotting/irregulardata3D.csv"),
source=csv,datatype=float[8]):
|
|
>
|
Ppt := pointplot3d(M,axes=box,style=patchnogrid,symbolsize=15,color=red,
labels=[x,y,z]):
|
|
>
|
display(Ppt,orientation=[90,0,0]);
|
ScatterPlot3D コマンドのデフォルトオプションでは、次数 2 の適合二次方程式を使用して、x 方向と y 方向の両方で元の領域の 1/3 幅のウィンドウを使用するように設定されています。
右下が平滑化した面、左が同じ面に元のデータポイントのポイントプロットを重ね合わせた透過表示です。
|
>
|
Ploess := ScatterPlot3D( M, lowess ):
|
|
>
|
display(Array([
display([Ppt,Ploess]),
display(Ploess,style=patch,color=RGB(0.0,0.4,0.6))
]),view=700..860,axes=box,orientation=[54,78,-15],
transparency=0.0);
|
|
|
|
|
補間
|
|
このセクションでは、1 次元の依存データが正確であることが判明している、2 次元独立データの集合があるものとします。目的は、すべてのデータポイントを通過する面を生成することです。つまり独立データの各 x-y ポイントで、プロット対象の面の高さが従属 (z) データの対応する値と一致しなければならないということです。
この場合、データを補完することで、面の数値計算が実行されます。
これは、2 つの異なる状況に分類されます。最初のケースでは、独立したすべての x-y データが規則的なグリッドを形成します。2 番目のケースでは、独立したデータポイントが x-y 平面上で不規則に配置されています。どちらのケースでも、surfdata コマンドを使って面をプロットできます。
|
|
均一なデータグリッド
|
|
均一な 2-D データグリッドの場合、面のプロットは、surfdata コマンドを使って生成できます。
データは x ポイントと y ポイントのグリッドから構成されます。このセクションの例では、データポイントを各方向で均一に取得します。x の値は1 ~ 7 の整数、y の値は 1 ~ 9 の整数として取得します (不均一なデータポイントについては、次のセクションを参照してください)。
|
|
x-y 平面のデータポイントのビュー
|
|
これは、入力データポイントのビューです。この例の場合、このビューは完全なグリッドで、x 方向と y 方向の両方で均一に配置されています。
|
>
|
xpts:=<($1..7)>:
ypts:=<($1..9)>:
|
|
>
|
display( seq(plottools:-line([xpts[i],1],[xpts[i],9]), i=1..7),
seq(plottools:-line([1,ypts[j]],[7,ypts[j]]), j=1..9),
pointplot([seq(seq([xpts[i],ypts[j]],j=1..9),i=1..7)], color=red) );
|
|
各データポイントに対応する指定したデータ値は、7 x 9 の配列のデータに保存されます。この例では、値の集合はランダムに生成されます。
|
>
|
data3D:=Array( LinearAlgebra:-RandomMatrix(7, 9, generator=0.2 .. 0.8), datatype=float[8]):
|
以下に示す最初の 3D プロットは、デフォルトで生成される区分関数 - 平面です。よりスムーズな面を生成することが目的で、これは、surfdata コマンドの gridsize と interpolation オプションを使用することで実現できます。
プロットされる面は、3D ポイントプロットと (線形の) パッチ面の両方と重ねて表示されます。このため、計算で得られた面が、元の離散データに対して確度が高いことがビジュアルに示されます。
|
>
|
ptsplot:=PLOT3D( GRID(1..7, 1..9, data3D, STYLE(POINT) ) ):
gridplot:=PLOT3D( GRID( 1..7, 1..9, data3D, STYLE(WIREFRAME)) ):
|
以下の plot3d オプション は、これらのプロットの一般的な外観をコントロールします。これらのオプションは、このサブセクションで再使用されます。 これらのオプションには、手法間の違いがより明確になるように、再利用可能な単一の名前が割り当てられています。
|
>
|
lookandfeel := axes=boxed, style=surface, glossiness=1.0, lightmodel=light1,
color=RGB(204/255, 84/255, 0/255), view=[1..7, 1..9, 0..1]:
|
|
>
|
display(
ptsplot,gridplot,
surfdata(data3D, 1..7, 1..9, lookandfeel)
);
|
左下は、surfdata コマンド (補間追加なし) によって生成されたデフォルトの面です。
右下は、新しい補間オプションであるキーワードパラメータ gridsize と interpolation を使った surfdata コマンドによって生成された面です。
gridsize パラメータは、gridsize=list という形式の方程式です。list には、元のデータの補間の実行対象となるポイントの均一グリッドのサイズを表す正の整数ペアを記述します。
interpolation パラメータは interpolation=list という形式の方程式です。list には CurveFitting:-ArrayInterpolation コマンドで使用できるオプションを記述します。
|
>
|
display( Array([
surfdata( data3D, 1..7, 1..9, lookandfeel ),
surfdata( data3D, 1..7, 1..9,
gridsize=[28,36], interpolation=[method=spline],
lookandfeel)
]) );
|
|
>
|
display( Array([
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9, lookandfeel ) ),
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9,
gridsize=[28,36], interpolation=[method=spline],
lookandfeel) )
]) );
|
3D プロットの表示は、少量単位による、面の独自平滑化によって実行されるという点に注意してください。したがって計算補間スキームという純粋に線形な手法を使用している場合でも、右下のプロットでは、適度なレベルの面の平滑化が実現されていることがわかります。
|
>
|
display( Array([
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9, lookandfeel ) ),
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9,
gridsize=[28,36], interpolation=[method=linear],
lookandfeel) )
]) );
|
値の補間を実行した精緻グリッドの場合、データポイントの元のグリッドの番号ポイントの全整数倍数を (両方向で) 含む必要はありません。
|
>
|
display( Array([
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9, lookandfeel ) ),
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9,
gridsize=[33,33], interpolation=[method=cubic],
lookandfeel) )
]) );
|
|
|
補間プロシージャ
|
|
次に、補間プロシージャ f3 と g3 を生成して使用します。各プロシージャは、指定 (x,y) ポイントに対して、補間処理したスカラー値を返します。このため、必要なポイント (x,y) を選択して、plot3d コマンドを直接使用することが可能になります。現在、plot3d は適応プロットを実行しないため、このような 2 パラメータ補間プロシージャを使用する利点はありません。
またこの手法は、(x,y) ポイントの事前作成済みグリッドを使用するという以前の方法と比べると、効率も低下します。これは ArrayInterpolation が起動され、入力ポイントごとに一次データ構造が生成されるためです。この件に関する効率の詳細については、付録を参照してください。
生成された補間プロシージャは、数値最適化、積分、根の計算など、他の一般計算にも使用できます。このような計算結果は、選択した補間手法によって異なるということを忘れないでください。
|
>
|
f3:=(x,y)->CurveFitting:-ArrayInterpolation([xpts,ypts], data3D,
Array(1..1, 1..1, 1..2, [[[x,y]]]), method=spline)[1,1]:
g3:=(x,y)->CurveFitting:-ArrayInterpolation([xpts,ypts], data3D,
Array(1..1, 1..1, 1..2, [[[x,y]]]), method=cubic)[1,1]:
drawfun3 := f_ ->
plot3d(f_, 1..7, 1..9, axes=boxed, style=surface,
glossiness=1.0, lightmodel=light1,
color=RGB(204/255, 84/255, 0/255)):
|
|
>
|
display(Array([drawfun3(f3), drawfun3(g3)]), view=[1..7, 1..9, 0..1]);
|
|
>
|
Optimization:-Minimize(f3, method=nonlinearsimplex);
|
|
![[-1.98108713564723402*10^61, Vector(2, {(1) = -0.9628329919e11, (2) = -5503890694.})]](/support/helpjp/helpview.aspx?si=2547/file01547/math682.png)
| (2.2.1.2.1) |
|
>
|
fsolve({(a,b)->a-b, (a,b)->f3(a,b)-0.5}, [1..7, 1..9]);
|
|
![[3.3483539130504871, 3.3483539130504871]](/support/helpjp/helpview.aspx?si=2547/file01547/math689.png)
| (2.2.1.2.2) |
|
| (2.2.1.2.3) |
|
|
|
|
不均一なデータグリッド
|
|
これは不規則な間隔で配置されたデータで、n x 3 の配列または行列として提供します。各行はポイントを表していて、列は x、y、または z の値と解釈されます。
次に例を示します。最初に、x-y 平面へのデータポイントの投影について説明するため、2D 従属データ値のレイアウトと間隔を表示します。
これは「平滑化」のセクションの ScatterPlot3D で使用したデータと同じです。ただし、平滑化の例と対比するため、データポイントの補間を実行するものとします。つまり、面がデータポイントを直接通過するということです。
|
>
|
M:=ImportMatrix(cat(kernelopts(mapledir),"/data/plotting/irregulardata3D.csv"),
source=csv,datatype=float[8]):
|
|
>
|
Ppt := pointplot3d(M,axes=box,style=patchnogrid,symbolsize=15,color=red,
labels=[x,y,z]):
|
|
>
|
display(Ppt,orientation=[90,0,0]);
|
surfdata コマンドは、このタスクを直接実行できるようになりました。
source=irregular オプションは、2 つの独立したデータ次元の 1 つに沿って 3 つの値しか持たない規則配置グリッドを形成するデータとして、このデータを解釈してはならないことを示すために指定します。
|
>
|
surfdata(M,source=irregular,axes=box);
|
|
|
|
|
|
付録 : 効率
|
|
CodeTools:-Usage コマンドは、個々の計算またはステップで使われるメモリリソースや時間リソースに関する情報を出力するために使用します。
効率に関する懸念の 1 つは、回避可能にもかかわらず、補間処理によって不必要なデータのコピーが実行されるということです。実際の作業を行うルーチンはすべて、プリコンパイルされた関数への外部呼出しとして実装されます。またこれらのルーチンでは、ハードウェア倍精度 datatype=float[8] オプションを指定した rtables (配列、ベクトル、または行列) として入出力を伝達することを前提としています。
1-D データに関する元の例を検討してみましょう。パフォーマンスの違いを明確に示すため、補間ポイント数 5000 という比較的大きな数を使用します。
出力ポイントの事前割当て構造 (A1) を使うと、補間プロットを最も効率的に取得できます。構造 A1 自体を生成するために必要なリソースは無視できます
(プロットの外観や平滑度については特に問題にしません。ここでは相対的な効率について調べます。個々にプロットされるポイントの数は必然的に大きくなり、場合によって二次的なジャギーが発生することがあります)。
|
>
|
A1:=CodeTools:-Usage(Vector(5000, i->i*15/5000, datatype=float[8])):
|
|
memory used=236.03KiB, alloc change=0 bytes, cpu time=0ns, real time=5.00ms
| |
|
>
|
CodeTools:-Usage(plot(<A1|ArrayInterpolation(data2D, A1, method=spline)>,
color=black));
|
|
memory used=27.77MiB, alloc change=0.59MiB, cpu time=340.00ms, real time=338.00ms
| |
F などの補間関数による、独立データの不必要な繰り返しコピーは、rtables を指定して ArrayInterpolation を実行することでも回避できます。
また ArrayInterpolation コマンドのオプションである container 引数を使用して、インプレースな動作を実現することで、スカラー出力結果を伝えるだけのために、rtable の作成とメモリ管理が繰り返し実行されるのを回避することができます。
最後に、補間関数をモジュールとして実装することで、従属入力ポイントを再利用可能な rtable として通知し、別の rtable オブジェクトのメモリ管理の繰り返しを回避することもできます。以下のモジュール Fnew は、このような概念を示しています。
|
>
|
Fnew := module()
export ModuleApply;
local inArray, outArray, hwxdata, hwydata;
inArray:=Array(1..1,'datatype'=float[8]);
outArray:=Array(1..1,'datatype'=float[8]);
hwxdata,hwydata := Vector(data2D[1..-1,1],'datatype'=float[8]),
Vector(data2D[1..-1,2],'datatype'=float[8]):
ModuleApply:=proc(x)
UseHardwareFloats:=true;
inArray[1]:=x;
CurveFitting:-ArrayInterpolation(hwxdata,hwydata,
inArray, 'method'='spline', 'container'=outArray);
outArray[1];
end proc;
end module:
|
ここで、 F と Fnew のパフォーマンスの違いを比較します。この例では、速度が約 2 倍改善します。スムーズで美しいプロットを実現するためには、5000 ものポイントは必要ないということを確認してください。 これほど多くのポイントを要求した目的は主に、パフォーマンスの違いを示すためです。
以下の 2 つのプロットは、出力ポイント数が同じで、事前作成済みの固定ベクトル A1 を使用した以前のプロットと比較して、かなり速度が低下しています。
|
>
|
CodeTools:-Usage( plot(F, 0..15, color=black, numpoints=5000) );
|
|
memory used=245.00MiB, alloc change=113.44MiB, cpu time=2.82s, real time=2.83s
| |
プロシージャ Fnew は F と比べて約 2 倍高速です。プログラミング工数を有意義に追加することで、Fnew による ArrayInterpolation の繰り返し呼び出しに起因するオーバヘッドを回避し、大幅な速度改善を実現することができます。F はスレッドセーフですが、Fnew はスレッドセーフではありません。これは、モジュール Fnew の inArray (または outArray) ローカル変数などモジュールローカル変数が、スレッドグローバルな手法でアクセスされるためです。Fnew に対する複数の個別コールが、異なるスレッドで同時に実行されると、こうしたローカル配列への同時アクセスは相互に干渉し合います。
|
|
