実際的な利用の例


ここでは実際的な利用の例として,ファイル(ファイル名はdata)から データを読み込み,処理を行う方法を紹介します.

データの読み込み

まず,ファイルdataの内容を確かめてみます. ファイルの内容は !!data で表示できます.

In[31]:=
!!data
大きなファイルの場合,次のようにすれば最初の10行だけを表示できます (Unixを利用している場合).

In[32]:=
ColumnForm[ReadList["!head -10 data", String]]
Out[32]=
0.1297323638545699
0.1323127172826521
0.1446228546294113
0.183239867688719
0.198061867847468
0.05509711267392398
0.1985483328594587
0.1346524383988073
0.1889659938776106
0.1923841700160055
一般に,ColumnForm[ReadList["!command", String]] によって, Unixコマンドの実行結果を表示できます. ReadListはファイルからデータを読み込むための関数ですが, ファイル名が ! から始まっている場合は,以降をUnixのコマンドと 解釈して,その実行結果を読み込みます(読み込んだ結果はリストとして 返されます). Stringは各行を文字列として読み込むための指定です. ColumnFormによりリストの要素を縦に並べて表示します.

Unixコマンドをしばしば実行する場合は, 次のような関数を定義しておくと便利です.

In[33]:=
Unix[command_String] :=
	Scan[Print, ReadList["!"<>command<>"|expand", String]]
In[34]:=
Unix["wc -w data"]
Out[34]=
     128 data
ファイルからデータを読み込むにはReadList関数を使います.

In[35]:=
data = ReadList["data", Number];
ReadListは,ファイルdataから数値データを読み込み, それをリストにして返しますので, そのリストを変数dataに代入しています.

各行に複数のデータがある場合の例です.次のようなファイルがあるとします.

In[36]:=
!!data2
Out[36]=
1	2	3
4	5
6	7	8
上と同じようにReadListを実行すると,単にすべての数値データを 一列にならべたリストが返ってきます. RecordListsオプションを指定すれば各行をリストにまとめることができます.

In[37]:=
data2 = ReadList["data2", Number]
Out[37]=
{1, 2, 3, 4, 5, 6, 7, 8}
In[38]:=
data2 = ReadList["data2", Number, RecordLists->True]
Out[38]=
{{1, 2, 3}, {4, 5}, {6, 7, 8}}
次のようにデータがコンマで区切られている場合はもう少しやっかいです. この場合,データの種別をWordと指定して,一旦データを文字列として 読み込んで,ToExpression関数により数値に変換する必要があります.

In[39]:=
!!data3
Out[39]=
1, 2, 3
4, 5
6, 7, 8
In[40]:=
data3 = ToExpression[
	ReadList["data3", Word, RecordLists->True,
			  WordSeparators->{" ", "\t", ","}]]
Out[40]=
{{1, 2, 3}, {4, 5}, {6, 7, 8}}
読み込んだデータを表示してみます.

In[41]:=
data
関数Shortをもちいれば,簡略化した表示をします.

In[42]:=
data // Short
Out[42]=
{0.1297323638545699, <<126>>, 0.1183108049642118}
expr // func は,func[expr] と全く同じことです. 上の場合は,Short[data] と同じことになります.

常に簡略化した表示にしたければ

In[43]:=
$PrePrint = Short
Out[43]=
Short
とします.次のようにすれば解除できます

In[44]:=
$PrePrint =.
読み込んだデータの平均値を計算してみましょう.

次のようにStatistics`DescriptiveStatistics`パッケージを 利用しても求められますが,直接求めるのも簡単にできます.

In[45]:=
Needs["Statistics`DescriptiveStatistics`"]
In[46]:=
Mean[data]					(* パッケージを利用 *)
Out[46]=
0.397853
In[47]:=
Plus @@ data / Length[data]		(* 直接求める *)
Out[47]=
0.397853
Plus @@ data は,Apply[Plus, data] と同じです. Apply[f, {x1,x2,...,xn}] は,f[x1,x2,...,xn] を求めます. つまり,今の例ではdataの要素の和が求められます. Plus[data] ではダメなのに注意して下さい. Length[data] は,リストdataの長さを返します. つまり,今の例ではdataの要素の個数が求められます.

データのプロット

では,読み込んだデータをプロットしてみます.

In[48]:=
ListPlot[data]
Out[48]=
-Graphics-
関数ListPlotは,リスト {y1, y2, ..., yn} が与えられたとき, すべての点 (i, yi) をプロットします.

実際の x-座標が1, 2, 3, ... なら良いのですが, そうでない場合はx座標とy座標の対のリスト {{x1,y1}, {x2,y2}, ..., {xn,yn}} を ListPlotに与えてやる必要があります. 今,xi = x0 + (i-1)dx だとすると,x-座標なしのリストから x-座標付きのリストへの変換は以下のようにしてできます. Table[expri, {i, n}] は, リスト {expr1, expr2, ..., exprn} を返す関数です. data[[i]] は,リストdataの第i要素を表します.

In[49]:=
Table[{x0+(i-1)dx, data[[i]]}, {i, Length[data]}]
しかし,次の方法のほうがエレガントです.

In[50]:=
Transpose[{x0+dx Range[0,Length[data]-1], data}]
まず,関数Rangeにより,リストdataと同じ長さの0から始まる リストを作成し,各要素をdx倍してx0を加えています. それとdataを並べたリストの転置行列が求める結果です.

これを関数として定義してみましょう.

In[51]:=
AddXcoord[list_, x0_, dx_] :=
	Transpose[{x0+dx Range[0,Length[list]-1], list}]
関数の定義は function_pattern := function_body のように, := を用います = は通常は用いません). また funciton_pattern 中の変数には必ず下線を付けてください.

:= と = の違いは,右辺を評価するかどうかの違いです. = は右辺を評価した結果を左辺に代入します. := は右辺を評価せずに左辺に代入します. 関数が実行されたときに初めて右辺を評価します. 以下の例を見て下さい.

In[52]:=
a = 10
Out[52]=
10
In[53]:=
f1[x_] := x+a	(* x+a が代入されます *)
In[54]:=
f2[x_] = x+a	(* x+10 が代入されます *)
Out[54]=
10 + x
In[55]:=
a = 20			(* aを20に変更 *)
Out[55]=
20
In[56]:=
f1[x]
Out[56]=
20 + x
In[57]:=
f2[x]
Out[57]=
10 + x
引数のデフォールト値を与えることもできます. AddXcoordで引数dxが与えられなければ1とするには, 以下のようにします.

In[58]:=
AddXcoord[list_, x0_, dx_:1] :=
	Transpose[{x0+dx Range[0,Length[list]-1], list}]
In[59]:=
AddXcoord[{10,11,12}, 1]
Out[59]=
{{1, 10}, {2, 11}, {3, 12}}
x-座標が -2 から 1/32 きざみだとして,実行してみます.

In[60]:=
datax = AddXcoord[data, -2, 1./32]
x-座標を取り除くのは簡単です.

In[61]:=
Transpose[datax][[2]]
では,x-座標をつけ加えたデータdataxをプロットしてみましょう.

In[62]:=
ListPlot[datax]
Out[62]=
-Graphics-
軸がグラフの真ん中にあるのが気になりますので, AxesOrigin->{-2,0} で軸の原点座標を指定します. ついでに AspectRatio->Automatic で,x-座標とy-座標のきざみを 同一にします(Mathematica は通常,グラフの外形を自動的に整形し, 黄金分割比を持つ横長の長方形になるようにします).

In[63]:=
ListPlot[datax, AxesOrigin->{-2,0}, 
				AspectRatio->Automatic]
Out[63]=
-Graphics-
図が小さい場合は,図をマウスでクリックして,図の枠をマウスでつかんで 引き伸ばしてください.

最小2乗法

Fit

データを最小2乗法で多項式近似するには,関数Fitを用います. ここでは,dataxを a + b x^2 + c x^4 + d x^6 という形の式に 当てはめ,求めた多項式を関数fit1と定義します.

In[64]:=
fit1[x_] = Fit[datax, {1,x^2,x^4,x^6}, x]
Out[64]=
                 2             4              6
1.034 - 1.39553 x  + 0.628783 x  - 0.0861296 x
一般の6次多項式で近似するには,Fitの第2引数を {1,x,x^2,x^3,x^4,x^5,x^6} または Table[x^i, {i,0,6}] としてください.

求めた多項式は Plot[fit1[x], {x,-2,2}] でプロットできますが, ここでは,もともとのデータdataxのプロットと重ねて表示してみます. 複数のプロットを重ねて表示するには, Graphics`Graphics`パッケージにあるDisplayTogether関数を使います.

In[65]:=
Needs["Graphics`Graphics`"]		(* パッケージの読み込み *)
In[66]:=
DisplayTogether[
		ListPlot[datax],
		Plot[fit1[x], {x,-2,2}],
		AxesOrigin->{-2,0}, AspectRatio->Automatic
]
Out[66]=
-Graphics-
次のようにしても複数のプロットを重ねて表示できます.

In[67]:=
graph1 = ListPlot[datax, DisplayFunction->Identity];
graph2 = Plot[fit1[x], {x,-2,2}, 
				DisplayFunction->Identity];
Show[graph1, graph2, AxesOrigin->{-2,0}, 
		AspectRatio->Automatic, 
		DisplayFunction:>$DisplayFunction]
Out[67]=
-Graphics-
DisplayFunction->Identityはプロットを表示しないための指定です. あとからShowのなかで DisplayFunction :> $DisplayFunction とすれば表示できます.

以降,たびたび重ねた表示を行うので,関数plotbothとして定義しておきます. plotboth[fit1]で上と同じプロットが得られます. ついでに2乗誤差を求める関数squareErrorも定義しておきましょう.

In[68]:=
plotboth[f_] := Module[{x},
	DisplayTogether[
		ListPlot[datax],
		Plot[f[x], {x,-2,2}],
		AxesOrigin->{-2,0}, AspectRatio->Automatic
	]]
squareError[f_] := Module[{i},
	Sum[N[(datax[[i,2]]-f[datax[[i,1]]])^2], 
		{i,1,Length[datax]}]
]
In[69]:=
squareError[fit1]
Out[69]=
0.743927
Module[{x1, x2, ..., xn}, body] は,局所変数x1, x2, ..., xnを宣言します (局所変数とはbodyの中だけで利用する変数です). plotbothでModuleを使わなかった場合,変数xはグローバル変数となり, xに値が代入されているときに正しく動きません.

NonlinearFit

多項式では,あまり良い近似が得られませんでした. 非線形の式で近似を行う場合にはStatistics`NonlinearFit` パッケージが使えます.

In[70]:=
Needs["Statistics`NonlinearFit`"]   (* パッケージの読み込み *)
In[71]:=
fit2a[x_] = a + b Exp[-(c x)^2];    (* あてはめる式を定義 *)
In[72]:=
NonlinearFit[datax, fit2a[x], x, {a,b,c}]
Out[72]=
{a -> 0.11288, b -> 1.01213, c -> 1.57377}
In[73]:=
fit2[x_] = fit2a[x] /. %            (* これが求める近似 *)
Out[73]=
            1.01213
0.11288 + -----------
                    2
           2.47676 x
          E
In[74]:=
plotboth[fit2]                      (* プロットは省略 *)
In[75]:=
squareError[fit2]
Out[75]=
0.422073
NonlinearFitは,Solveなどと同様に解を規則の形で返しますので,
/. を使って規則を適用した結果を求める必要があります.

NonlinearFitを使わなくても,y-座標データの Logを 計算すれば,Fitでももう少し良く近似できます.

In[76]:=
dataxLog = Map[{#[[1]], Log[#[[2]]]}&, datax];
In[77]:=
fit3[x_] = Exp[Fit[dataxLog, {1,x,x^2,x^3,x^4}, x]]
Out[77]=
                                           2
Power[E, 0.0576518 - 0.184933 x - 2.06158 x  + 
 
              3            4
   0.0110031 x  + 0.40145 x ]
In[78]:=
plotboth[fit3]                      (* プロットは省略 *)
In[79]:=
squareError[fit3]
Out[79]=
0.5461

FFT

最後に離散型フーリエ変換で近似してみましょう. まず次のような関数を定義します(説明は省略).

In[80]:=
FitFFT[data_, {x_,x0_:0,dx_:1}, m_] := Module[{n,f,k},
	n = Length[data];
	f = N[InverseFourier[data] / Sqrt[n]];
	Re[f[[1]]] +
	   Sum[ 2 Re[f[[k+1]]] Cos[2 Pi k (x-x0) / (n dx)] -
		    2 Im[f[[k+1]]] Sin[2 Pi k (x-x0) / (n dx)] ,
            {k,1,m} ]
]
3次まで近似した式を求めます.

In[81]:=
fit4[x_] = FitFFT[data, {x,-2,1/32}, 3]
Out[81]=
                        Pi (2 + x)
0.397853 - 0.442581 Cos[----------] + 
                            2
 
  0.216662 Cos[Pi (2 + x)] - 
 
                3 Pi (2 + x)
  0.0506513 Cos[------------] + 
                     2
 
                Pi (2 + x)
  0.0272171 Sin[----------] + 
                    2
 
  0.0056861 Sin[Pi (2 + x)] + 
 
                 3 Pi (2 + x)
  0.00208212 Sin[------------]
                      2
In[82]:=
plotboth[fit4]
Out[82]=
-Graphics-
In[83]:=
squareError[fit4]
Out[83]=
0.378544

データの書き出し

最後に,データをファイルに書き出すプログラムを紹介します. Mathematica だけで利用するデータの場合は,以下のようにして 簡単にファイルに保存できます.

In[84]:=
datax >> "tmp1"           (* ファイル tmp1 に保存 *)
In[85]:=
newdatax = << "tmp1"      (* ファイル tmp1 から読み込み *)
In[86]:=
datax == newdatax	(* dataxとnewdataxは等しい *)
Out[86]=
True
ただし,tmp1には Mathematica の形式で書き込まれるので, Mathematica 以外のプログラムで取り扱うのは面倒です.

In[87]:=
!!tmp1
このような時は,次のプログラムを利用すればうまくいくと思います (説明は省略).

In[88]:=
WriteArray[file_, array_List] := Module[{stream},
	stream = OpenWrite[file, PageWidth->Infinity];
	Scan[WriteToStream[stream, #]&, array];
	Close[stream];
]
WriteToStream[stream_, list_List] := Block[{},
	If[ Length[list] > 0,
		WriteString[stream, First[list]];
		Scan[WriteString[stream, "\t", #]&, Rest[list]]
	];
	WriteString[stream, "\n"];
]
WriteToStream[stream_, x_] := Write[stream, x]
たとえば,dataxは以下のような形式で保存されます.

In[89]:=
WriteArray["tmp2", datax]
In[90]:=
!!tmp2

最初に戻る 前へ
Dept. CS / Faculty of Eng. / Kobe Univ. / Naoyuki Tamura / tamura@kobe-u.ac.jp