統計計算電卓

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() 関数で計算します。次の例は3列の乱数のデータを発生させて、各列の間の相関係数を計算しています。

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 のように関数名を指定すると使い方の簡単な説明を見ることができます。