無相関検定のシミュレーション

2種類のデータの間に相関が無いとき、標本相関係数rを使って定義した
\[ t = \sqrt(n-2)\frac{r}{\sqrt{1-r^{2}}}\]
は自由度n-2のt分布に従うことが、理論的にわかります。各データが母集団から無作為抽出されたとき、正規分布に従って標本が抽出されるとすれば、この量がt分布に近づくことをR言語を使って次のように確かめます。

#test_cor.R
size <- 20000
sample_size <- 20
mu1 <- 50
sigma1 <- 10
mu2 <- 30
sigma2 <- 20
ts <- numeric(length=size)
for(i in 1:size) {
  sample_x <- rnorm(n=sample_size, mean=mu1, sd=sigma1)
  sample_y <- rnorm(n=sample_size, mean=mu2, sd=sigma2)
  xbar <- mean(sample_x)
  ybar <- mean(sample_y)
  sxy <- mean((sample_x - xbar)*(sample_y - ybar))/sample_size
  sxx <- mean((sample_x - xbar)^2)/sample_size
  syy <- mean((sample_y - ybar)^2)/sample_size
  r <- sxy/sqrt(sxx*syy)
  ts[i] <- sqrt(n-2)*r/sqrt(1-r^2)
}

Rのコンソールから次のようにして、図を書き出します。

> png("test_cor.png")
> hist(ts, freq=FALSE)
> curve(dt(x, 18), add=TRUE)
> dev.off()

f:id:nabeyang:20130328090339p:plain

母分散が未知のときの平均値の検定のシミュレーション

母分散が未知のときの平均値の検定に使う検定量はサンプルの大きさをnとすると、自由度n-1のt分布に従うことが計算から分かります。これをR言語で、以下のコードで確かめることができます。このコードではn=20の場合にあたります。

# test_t.R
size <- 10000
sample_size <- 20
mu <- 50
sigma <- 10
ts <- numeric(length=size)
for(i in 1:size) {
  sample <- rnorm(n=sample_size, mean=mu, sd=sigma)
  xbar <- mean(sample)
  z <- (xbar - mu)/sqrt(sigma^2 /sample_size)
  y <- sum((sample - xbar)^2/sigma^2)
  ts[i] <- z/sqrt(y/(sample_size-1))
}

これを、自由度19のt分布と比較します。

> png("sample_t.png")
> hist(ts, freq=FALSE)
> curve(dt(x, 19), add=TRUE)
> dev.off()

f:id:nabeyang:20130327175808p:plain

不偏分散のバラつきのシミュレーション

標本の大きさnの不偏分散のバラつきは、(n-1)(不偏分散)/(母分散)が理論的に自由度(n-1)のカイ二乗分布に従うことから分かります。ここではR言語で、そのことをシミュレーションで確かめる方法を紹介します。
まずtest_chi2.Rというスクリプトを書きます。これはn=10の(n-1)(不偏分散)/(母分散)を10000回計算した結果をvarsに記録しています。

#test_chi2.R
size = 10000
vars <- numeric(length=size)
sigma <- 10
for(i in 1:size) {
  sample <- rnorm(n=10, mean=50, sd=sigma)
  mu <- mean(sample)
  vars[i] <- sum(((sample - mu)/sigma)^2)
}

これをRのコンソールから呼び出して、グラフを書かせます。

> png("sample_chisq.png")
> hist(vars, probability=TRUE)
> curve(dchisq(x, 9), add=TRUE)
> dev.off()

結果は次のようになります。
f:id:nabeyang:20130327101814p:plain

ソートアルゴリズムを暗記するためのゲームを作ってみた

最近、courseraでアルゴリズムを勉強していました。そのときやった練習問題の1つのタイプに「xxxソートで4回スワップした結果を書きなさい」というものがあって、これがなかなか勉強になりました。イメージができれば、コーディングはだいたいできるので、このような作業を繰り返しやれば、空で基本的なソートアルゴリズムくらいコーディングできるようになると思います。

それで、この練習問題をするとき、紙に書いてやっていくとなかなか時間がかかります。数回の操作につき1回は全ての配列の状態を書き直す必要があるからです。それで、もうちょっと気楽にポチポチとやれるものが欲しいなと思ってゲームっぽいものを作ることにしました。

それでできたのが、これです(ただしまだコマンドライン上でしか遊べない)。

ゲームっぽいところは、スコアがあるところと、ミスせずに成功した回数がコンボとして数えられるところ、コンボが続くとボーナスポイントが入るところです。普通に起動すると、Selection, Insertion, Shell, Heap, Quick, Mergeでランダムに与えられた長さ10の配列をソートすることになります。ソートし終わると、新しい配列が生成されます。

作った動機

  • 紙で書くのは大変。カードでやっても良いけれどもどれくらい頑張ったか見れない
  • なんらかの試験を受けるときは、暗記しないと話にならない
  • なんらかのアルゴリズムを設計するとき、基本的なアルゴリズムを理解していないと話にならない
  • ゲームとしても面白いんじゃないか?
    • いくつかのタイプの作業を言われたり、見て瞬間的に判断してやることを繰り返すようなタイプのゲームはけっこうある

それぞれのモードの動画

Selection

mrkrで"1"と表示されているindexのところまでソートされています

Insertion

Shell

間隔はh = 3 * h + 1で決めます

Heap

配列をツリー状にしたものも表示します

Merge

Quick

matplotlibを使って計算時間をプロットしてみた

matplotlibを使って、JavaのHashSetとArrayListのcontainsの実行速度を測った結果を図にしました。loglogグラフなので、ArrayListはO(N)でHashSetはO(1)でしょうか。測定方法とテストの詳細は「実装パターン」の付録Aに載っています。

使ったデータ

N 1 10 100 1000 10000 100000
setMembership 7.26 9.97 10.49 10.94 11.52 11.21
arrayListMembership 6.65 29.72 228.36 2585.26 35996.90 1257285.81

matplotlibのコード

#!/usr/bin/env python2.7
#encoding:utf-8
from matplotlib import pyplot
with open("set_arraylist.data") as f:
  x = [int(d) for d in f.readline().rstrip().split(" ")[1:]]
  ys = []
  legends = []
  for line in f:
    data = line.rstrip().split(" ")
    legends.append(data[0])
    ys.append([float(fl) for fl in data[1:]])
  size = len(x)
  for y in ys:
    if len(y) == size:
     pyplot.plot(x, y)
  pyplot.legend(legends)
  pyplot.xlabel("#elements")
  pyplot.ylabel("processing time(ns)")
  pyplot.xscale("log")
  pyplot.yscale("log")
  pyplot.show()

matplotlibのインストール

環境はUbuntu 12.04です。pipでも入りましたが、showでGUIが立ち上がりませんでした。色々試しましたが結局apt-getで一発で入りました。公式サイトに書いてあるのと同じ方法です。

$ apt-get install python-matplotlib

Javaオブジェクトの使用メモリの計算するプログラムを書きました

Javaの定義されたクラスから生成されたオブジェクトのメモリ使用量(64bitマシン)の計算をするプログラムを書きました。プリミティブなフィールドのみに対応しています。(オーバーヘッドの部分と参照のところが32bitとずれます。詳しくは、例えばここを参照してください。)。"gen_nomethod_class.py"はテスト用のクラスの生成のために使っています。

使い方は次のようになります。

$ java MemoryCheck  ClassName
usage = 4296 bytes

Javaで"Hello World"

Javaは(windowsで昔)Eclipseを使って書いたことはあったのですが、コマンドラインから実行したことはなかったので、調べました。環境はUbuntu 12.04で、ここに書いてある通りにしてjavaを導入しました。次に"Hello World"を書いて、実行方法を確認しました。

// HelloWorld.java
public class HelloWorld {
  public static void main(String [] args) {
    System.out.println("Hello, World!");
  }
}

を作って、

$ javac HelloWorld.java

とすると、HelloWorld.classが生成されるので

$ java HelloWorld

とすると"Hello, World!"と標準出力されます。