円周率1,000,000桁で各数字が等確率で出現するか問題

  • 0から9の各数字が、円周率1000桁に等確率で出現するかを仮説検定してみる

一般に、

ある属性Aによって、n個の個体がk種のカテゴリー A_1, A_2, ..., A_k へ分類されるとき
各カテゴリーへ属する観測度数が f_1, f_2, ..., f_k であるとする
これが、各カテゴリーの理論確率 p_1, p_2, ..., p_k に適合するかを見るには、これが正しければ生じるであろ
う理論度数 np_1, np_2, ..., np_k を、観測度数と比べ、K.ピアソンの適合度基準
\chi^2 = \Sigma^k_{i=1} \frac{(f_i-np_i)^2}{np_i}
で判断すればよい。
この適合度のχ^2統計量χ^2はnが大きいとき、自由度 k-1 の χ^2 分布 \chi^2(k-1) に従う。

帰無仮説
H_0: P(A_1)=p_1, ..., P(A_k)=p_k
とするとき、
\chi^2 > \chi^2_\alpha(k-1)
ならば、"観測度数は理論確率分布 p_1, ..., p_k に適合している" という仮説H_0は、有意水準\alphaで棄却される。ここで、自由度=カテゴリー数(k)-1である。

円周率1,000桁でやってみる

計算

# observed frequency of 0, 1, ..., 9 (これは適当にスクリプトでカウントした)
freq <- c(93,116,103,102,93,97,94,95,101,106)

# expected probability
p <- 0.1

# number of observations
n <- 1000 

np <- n * p

chi <- sum((freq - np)^2 / np)
# => 4.74 

qchisq(0.95, 9)
# => 16.91898

chiの値(4.74) < \chi^2(9)(16.91898) なので、有意水準5%で棄却されない(0, ..., 9 は等確率で出現している)。

1,000,000桁でもやってみる

# これは適当にスクリプトでカウントした
freq <- c(99959,99758,100026,100229,100230,100359,99548,99800,99985,100106)
p <- 0.1
n <- 10^6
np <- n * p

chi <- sum((freq - np)^2 / np)
# => 5.50908

qchisq(0.95, 9)
# => 16.91898

\chi^2(5.50908) < \chi^2_{0.05}(9)(16.91898) なので、棄却されない。

同じ事だけどp値の計算でやると

pchisq(5.50908, 9, lower.tail=F)
# => 0.7878669

p値(0.7878669) > 0.05 (有意水準) となるので帰無仮説は棄却されない。

リハビリ

「Web+DB vol.66 きっちりわかるアルゴリズム」の「埋め込み性的辞書による圧縮アルゴリズム」をHaskell

module Main where

type Dict = [(String, String)]

dict :: Dict
dict = [ ("http", "'")
       , (".jp", "|")
       , ("www", ">")
       , (".co", "<")
       , ("html", "^")
       ]

replace :: Eq a => [a] -> [a] -> [a] -> [a]
replace [] _ _ = []
replace s find repl =
    if take (length find) s == find
        then repl ++ (replace (drop (length find) s) find repl)
        else [head s] ++ (replace (tail s) find repl)

compressWithDict :: Dict -> String -> String
compressWithDict [] s = s
compressWithDict ((find, repl):rest) s = compressWithDict rest $ replace s find repl

decompressWithDict :: Dict -> String -> String
decompressWithDict d = compressWithDict dict'
    where dict' = [ (y,x) | (x,y) <- d ]

compress :: String -> String
compress = compressWithDict dict

decompress :: String ->String
decompress = decompressWithDict dict

統計勉強する

今日から統計とかRとかやる
ネタは

Rによるやさしい統計学

Rによるやさしい統計学

はじめての統計学

はじめての統計学

dfltweb1.onamae.com – このドメインはお名前.comで取得されています。
実践! Rで学ぶ統計解析の基礎 − @IT
Statistical Aspects of Data Mining (Stats 202) Day 1 - YouTube

AnyEventの練習

噂のsleepsort書いてみた。AnyEvent->timer の戻り値捨てると動かないんだなぁとか。

#!/usr/bin/env perl

use strict;
use warnings;
use Data::Dumper;
use AnyEvent;

sub sleepsort {
    my (@args) = @_;

    my $cv = AnyEvent->condvar;
    my (@sorted, @timers);

    for my $arg (@args) {
        $cv->begin;
        push @timers, AnyEvent->timer(
            after => $arg / 10,
            cb => sub {
                push @sorted, $arg;
                $cv->end;
            },
        );
    }

    $cv->recv;
    \@sorted;
}

print Dumper sleepsort(@ARGV);

Perlのsystem()に@_を渡せなくてハマった。と思ったら単純ミス。

背景

コマンドラインで動くPerlのプログラムをTDDで作ろうと思い、同じくコマンドラインプログラムであるApp::cpanminusのテストの書き方を真似して書いていた。

問題

でもユーティリティ関数で呼んでいるsystem()に引数が渡らない!なぜ!><

純化するとこんなプログラム

#!/usr/bin/env perl

use strict;
use warnings;
use Capture::Tiny qw(capture);

sub caller0 {
    my $out = capture { system('echo', @_) };
    print $out;
}

sub caller1 {
    my @args = @_;
    my $out = capture { system('echo', @args) };
    print $out;
}

print "caller 0\n";
caller0('output!');  # 何も出ない

print "caller 1\n";
caller1('output!');  # ちゃんと出力される

caller0内で直接@_を使うと動かない。引数に何も渡ってこない。

解決編

ひたすらハマったけど気づいてみると当たり前だった。captureは単なる関数なので、{ ... } は無名サブルーチンを渡しているだけ。なので caller0の引数のつもりで@_と書いても、実際はcaptureに渡しているサブルーチンの引数と解釈される。だからcapture呼ぶ前にコピーしてそっちを渡しましょうっていうだけの話。

おまけ

ちなみに参考にしてるApp::cpanminusのコードはこちら

xt/Run.pm

package xt::Run;
use strict;
use base qw(Exporter);
our @EXPORT = qw(run last_build_log);

use Capture::Tiny qw(capture);
use File::Temp qw(tempdir);

$ENV{PERL_CPANM_HOME} = tempdir(CLEANUP => 1);

sub run {
    my @args = @_;
    my @notest = $ENV{NOTEST} ? ("--notest") : ();
    return capture { system($^X, "./script/cpanm.PL", @notest, "--quiet", "--reinstall", @args) };
}

sub last_build_log {
    open my $log, "<", "$ENV{PERL_CPANM_HOME}/latest-build/build.log";
    join '', <$log>;
}

1;

最初これをコピペしてたのに何で書き換えちゃったんだろう。。。まぁ勉強になりました。

memcached

You can enable memcached to automatically load on login with:
    mkdir -p ~/Library/LaunchAgents
    cp /usr/local/Cellar/memcached/1.4.5/com.danga.memcached.plist ~/Library/LaunchAgents/
    launchctl load -w ~/Library/LaunchAgents/com.danga.memcached.plist

    Or start it manually:
        /usr/local/bin/memcached

    Add "-d" to start it as a daemon.