How to use Mathematica 6 (in Japanese)

実際的な利用の例

ここでは実際的な利用の例として，ファイル(ファイル名は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]=

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]=

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