Octave は統計計算も得意です。Octave を統計計算電卓として使ってみましょう。
floor() 関数(小数点下切捨て) と rand() 関数(乱数発生)を使ってサンプルデータを作ります。
$ octave GNU Octave, version 2.0.16 (i386-vine-linux-gnu). Copyright (C) 1996, 1997, 1998, 1999, 2000 John W. Eaton. This is free software with ABSOLUTELY NO WARRANTY. For details, type `warranty'. octave:1> a = floor( rand(9,1) * 100 ) a = 55 27 48 44 41 79 18 50 42
基本統計量は statistics() で計算します。
octave:2> statistics( a )
ans =
18.00000 ...... min (最小値)
41.00000 ...... first quartile (第一四分位)
44.00000 ...... median (中央値)
50.00000 ...... third quartile (第三四分位)
79.00000 ...... max (最大値)
44.88889 ...... mean (平均値)
17.20788 ...... std (標準偏差)
0.34102 ...... skewness (歪度)
-0.47753 ...... kurtosis (尖度)
min, median, max, mean, std, skewness, kurtosis などは単独でも計算できます。
octave:3> mean( a ) ans = 44.889 octave:4> std( a ) ans = 17.208
データの並べ換えは sort() 関数でできます。
octave:5> sort( a ) ans = 18 27 41 42 44 48 50 55 79
相関係数は corrcoef() 関数で計算します。次の例は３列の乱数のデータを発生させて、各列の間の相関係数を計算しています。
octave:6> b = rand(10, 3) b = 0.4790020 0.4872454 0.8961018 0.9393837 0.7651517 0.0738843 0.4490651 0.1958129 0.1467713 0.4672043 0.6830208 0.5775421 0.6191948 0.9705683 0.7874286 0.3321203 0.7263469 0.8015056 0.0097532 0.9273366 0.8649854 0.0319095 0.6078545 0.8053926 0.2698602 0.6702657 0.2115912 0.2568175 0.1534261 0.2287156 octave:7> corrcoef( b ) ans = 1.000000 0.096873 -0.411244 0.096873 1.000000 0.478464 -0.411244 0.478464 1.000000
分割表のχ二乗検定は chisquare_test_independence () 関数を用いて簡単に行えます。次の 2 × 2 分割表のχ二乗検定をしてみましょう。
octave:1> A = [10, 3; 5, 7] A = 10 3 5 7 octave:2> [pval, chisq, df] = chisquare_test_independence ( A ) pval = 0.072220 chisq = 3.2318 df = 1
この例ではχ^2 = 3.2318、自由度 1 で帰無仮説の棄却確率が 7.2 % なので分割表 A には有意な傾向性はないと結論づけられます。また、左辺が省略されると、棄却率の表示だけとなります。
octave:3> chisquare_test_independence ( A ) pval: 0.0722196 ans = 0.072220
二組のデータの平均値が有意に差があるかどうかを検定するのには t 検定を用います。
先ず最初に標準偏差が 1 で平均値がそれぞれ 0 と 1 からなる 20 個ずつのデータ x と y を次のように作ります。
octave:1> x = randn(20,1); octave:2> y = randn(20,1) + 1;
x と y の標本平均と標準偏差の推定値は次のようになります。
octave:3> mean( x ) ans = 0.60581 octave:4> std( x ) ans = 0.78542 octave:5> mean( y ) ans = 1.3597 octave:6> std( y ) ans = 1.0366
また、次のようにしてサンプルデータをグラフにプロットします。
octave:7> u = ones(20,1); octave:8> v = u + 1; octave:9> gplot [0:3][:] ([u,x]) with points, ([v, y]) with points
そこで、x と y の平均値に有意差があるかどうか t 検定します。
octave:10> [p, t, df] = t_test_2 (x, y, "<>") p = 0.013456 t = -2.5923 df = 38
両側検定では有意水準 1.3% で有意差があることが分かります。x の平均値が y の平均値より小さいと考えた場合の片側検定では次のようになります。
octave:11> [p, t, df] = t_test_2 (x, y, "<") p = 0.0067282 t = -2.5923 df = 38
この場合は、有意水準 0.6% で有意差があると判定されます。逆に x の方が大きいと考えた場合はどうでしょう。
octave:12> [p, t, df] = t_test_2 (x, y, ">") p = 0.99327 t = -2.5923 df = 38
この場合は、有意水準 99% で仮説は棄却されます。
正規分布の確率密度関数は normal_pdf (X, mean, std) で与えられます。次のように操作すると正規分布の確率密度関数をグラフで表すことができます。
octave:1> X = linspace(-1, 1, 100)' octave:2> Y = normal_pdf ( X, 0, 0.1) octave:3> plot (X, Y)
その外の確率分布は octave のプロンプト上で help と入力して help 画面を出したた後、/statistics で検索すれば捜すことができます。また、octave のプロンプト上で help normal_pdf のように関数名を指定すると使い方の簡単な説明を見ることができます。