Rule of Three

調べ物をしていて、このような法則があることをふと思い出した。

病害虫の発生率がp%以下であるということを95%の信頼度で言うためには、3/p個の対象を調べて病害虫が存在しないことを示せば良い。

この法則はRule of Threeあるいは3の法則と呼ばれている*1

Rule of Threeは医薬品の稀な副作用を調査するという文脈の上で、次のように言われることもある。おそらくこちらの表現の方が有名だろう。

100人調べて1度も観測されなくても別の100人中3人に事象が生起する可能性があり、100 万人調べて0でもなお他の100万人中3人に生起する可能性は否定できない

導出

なんとなく「ポアソン分布から導出できる」みたいなアバウトな覚え方をしていて、いざやってみたら意外と手こずってしまったのでメモしておく。

まず前提条件を整理すると次のようになる。

  • 事象が起こる確率pはすごく小さい(と予想している)
  • 試行回数nは大きい
  • n回の調査で1回も事象は発生しなかった

pは二項分布に従うと考えられるが、pが大きくnが小さい二項分布はnp=\lambdaをパラメータとするポアソン分布で近似できる。ポアソン分布において事象がk回起こる確率は

P(k) = \frac{\lambda^k e^{-\lambda}}{k!}

だが、今回k=0なので、

P(0) = e^{-\lambda}

となる。ここで、\lambdaの上限を考えてみよう。\lambdaが大きくなるほどPは小さくなるだけで\lambdaはいくらでも大きくできることはできる。ただ、あまり小さなPを想定するのは不自然なので、P=0.05に対応する\lambdaあたりを上限としよう。

わざとらしい感じになってしまったが要するに\lambdaの片側95%信頼区間を求めるということだ。

0.05 = e^{-\lambda}

\lambdaについて計算すると、

\lambda = -\ln{0.05} \approx 3

となる。要するに、n回の試行で事象の発生回数が0だったとしても、\lambda=n回の試行における事象の発生回数の期待値の95%信頼区間はほぼ0〜3ということになる。

これが

100 万人調べて0でもなお他の100万人中3人に生起する可能性は否定できない

に対応する。ここでさらにnp=\lambdaとしてnについて整理すると

n = \frac{3}{p}

となり、これが

病害虫の発生率がp%以下であるということを95%の信頼度で言うためには、3/p個の対象を調べて病害虫が存在しないことを示せば良い。

に対応する。

参考

*1:私はこの法則が植物の病害虫調査で実際に利用されている例は見たことがない。住んでいた分野が違うだけかもしれないが

連続した空行(空白文字を含む)を置換する正規表現

連続する空行の置換

O'Reilly Japan - 詳説 正規表現 第3版に、連続した空行(これは、任意の空白文字を含む場合がある)を<p>に置換する正規表現として、

$text =~ s/^\s*$/<p>/mg;

というものが掲載されている(p.67)。これは、空白文字を含む連続した空行を含む文字列...with.\n\n\n\t \nTherefore......with.\n<p>\nTherefore...のように置換すると説明されていて、実際この例ならうまくいく。

3つ以上連続する改行が正しく置換できない

しかし、a\n\n\nbのように改行が3つ以上続く例はa\n<p><p>\nbと置換される。改行がいくつでも<p>は2つだ。『詳説 正規表現』は10年近く前の本なので、Perlの仕様変更があったのかもしれないけど、ちょっと調べた程度ではどうしてこのような挙動になるのかはよくわからなかった。ちなみにa\n\n\n\t\n\nba\n<p><p>\nbに置換されたりしてますますよく分からない。

パターン修飾子\mを付けた場合に連続する改行がどう解釈されるのか曖昧とかそんなことが理由じゃないかと思っていたけど、どうやらこれは/gを付けてグローバルマッチをさせた場合に起こるらしかった。/gを付けた場合に実際に何が起こるか、という点を正確に理解できていないのかもしれない。

/gを指定せずs/^\s*$/<p>/mとすれば予想通りの動きをする(もちろん1箇所しか置換できない)。

ということは、

while($text =~ s/^\s*$/<p>/m){}

とかやってやれば期待通りの置換ができる。

別解

とはいえ/gがあったりなかったりで挙動が変わるものを使うのは気持ちが悪いので、他の方法を探していたところStack Overflowに答えを見つけた(cf. Perl Regex To Condense Multiple Line Breaks - Stack Overflow)。

Perl 5.10から導入された文字クラスとして、\R\hがある。\Rは改行っぽいものにマッチする文字クラス(cf.perlrebackslash - Perl 正規表現逆スラッシュシーケンスとエスケープ - perldoc.jp)で、\hは水平方向の空白文字(タブとかスペースとか)にマッチする文字クラス(cf. perlrecharclass - Perl 正規表現文字クラス - perldoc.jp)であって、これを組み合わせると、連続した空行の置換は、

$text =~ s/\R(\h*\R)+/\n<p>\n/mg;

と書ける。\R\hは使えない言語もあるかもしれないが、こちらの方が何をやっているのか明確という印象はある。

他の言語

ところで、他の言語でも同じような現象が起こるかというとそうではないっぽい。途中で面倒になったのであんまり試してないけど。。

R

例えばR、gsub()で試してみると問題なく置換できる。

Python3

PythonもOK。

Go

Goもいけた。

JavaScript

試した範囲でPerlと同じ現象が起こるのはJavaScriptだった。/gをオプションで指定するタイプだと発生するっぽい?

Project Euler Problem 24 - 順列の辞書順ソート

Problem 24 - Project Euler

順列とは要素の並び替えのことである。たとえば3124は数字1, 2, 3, 4の並び替えで得られる順列の一つだ。全ての順列をアルファベット順またはABC順に並び替えたとして、それを辞書順と呼ぶことにする。0, 1, 2の順列を辞書順に並び替えたとすると、結果は

012 021 102 120 201 210

となる。0, 1, 2, 3, 4, 5, 6, 7, 8, 9の順列を辞書順に並び替えたときの100万番目は何か。

以下解答。

続きを読む

Project Euler Problem 23 - 2つの過剰数の和ではない自然数の和

自分自身を除く約数の和が自分自身に等しい自然数完全数と呼ばれる(e.g. 6 = 1 + 2 + 3)が、そうでないものについて、不足数過剰数という概念がある。

  • 不足数: 自分自身を除く約数の和が自分自身より小さい自然数
  • 過剰数: 自分自身を除く約数の和が自分自身より大きい自然数

ある自然数を2つの過剰数の和で表すことを考える。例えば12は最小の過剰数(1 + 2 + 3 + 4 + 6 = 16)なので、2つの過剰数で表せる最小の自然数は24だ。

ところで、このような自然数は無限にある。それどころか、28123より大きな自然数全て2つの過剰数の和として表すことができることが知られている。

では、過剰数の和として表すことができない自然数の総和を求めよ、というのがこの問題。

最初bruteforceでやって約80秒だった。forum覗いて少し良いやり方を見つけたものの、劇的な改善はせず38秒。そしてやり方はわかるもののなぜこれで正しい答えが出るのかも理解しきれていない。計算時間については、他の言語の様子を見ると普通にネストしたfor回して200msとかになってるし、まあ仕方ないかなという気持ちもある。

ところで、28123という上限は20161まで下げられる。Wikipediaにも書いてある(過剰数 - Wikipedia)。ループがネストしているので計算量は2/3より大きく減らせて、計算時間は20秒になった。

以下解答。

続きを読む

Pythonによるハンバーガー統計学 #2

引き続きハンバーガー統計学にようこそ!の第2章。

準備

import numpy as np
import pandas as pd
import scipy as sp

2.1 平均的ポテトを推定する

wakupote = pd.Series([47, 51, 49, 50, 49, 46, 51, 48, 52, 49])
wakupote.describe()

2.3 信頼区間/区間推定

p値からt値への関数はscipy.stats.t.ppf

from scipy.stats import t
# 95%信頼区間に対応するt値
t.ppf(q=[0.025, 0.975], df=9)
# 同99%
t.ppf(q=[0.005, 0.995], df=9)
array([-2.26215716,  2.26215716])
array([-3.24983554,  3.24983554])

信頼区間scopy.stats.t.intervalで求める。dfに自由度、locに標本平均、scaleに標本標準誤差を指定する。標本標準誤差はpandas.Series.semで求められる(cf. Pythonで正規分布の平均値の信頼区間を計算する方法)。

t.interval(alpha=0.95, df=9, loc=wakupote.mean(), scale=wakupote.sem())
t.interval(alpha=0.99, df=9, loc=wakupote.mean(), scale=wakupote.sem())
(47.859567155669005, 50.540432844331001)
(47.274321990699256, 51.125678009300749)

自由度を複数与えることもできるので、t分布表も作れる。

df = np.hstack((np.arange(1,31,1), [40, 60, 120, 1e10]))
pd.DataFrame(data = {"p=0.95" : t.ppf(q=[0.975], df=df),
                     "p=0.99" : t.ppf(q=[0.995], df=df)},
             index=df)
                p=0.95   p=0.99
1.000000e+00    12.706205   63.656741
2.000000e+00    4.302653    9.924843
3.000000e+00    3.182446    5.840909
4.000000e+00    2.776445    4.604095
5.000000e+00    2.570582    4.032143
6.000000e+00    2.446912    3.707428
7.000000e+00    2.364624    3.499483
8.000000e+00    2.306004    3.355387
9.000000e+00    2.262157    3.249836
1.000000e+01    2.228139    3.169273
1.100000e+01    2.200985    3.105807
1.200000e+01    2.178813    3.054540
1.300000e+01    2.160369    3.012276
1.400000e+01    2.144787    2.976843
1.500000e+01    2.131450    2.946713
1.600000e+01    2.119905    2.920782
1.700000e+01    2.109816    2.898231
1.800000e+01    2.100922    2.878440
1.900000e+01    2.093024    2.860935
2.000000e+01    2.085963    2.845340
2.100000e+01    2.079614    2.831360
2.200000e+01    2.073873    2.818756
2.300000e+01    2.068658    2.807336
2.400000e+01    2.063899    2.796940
2.500000e+01    2.059539    2.787436
2.600000e+01    2.055529    2.778715
2.700000e+01    2.051831    2.770683
2.800000e+01    2.048407    2.763262
2.900000e+01    2.045230    2.756386
3.000000e+01    2.042272    2.749996
4.000000e+01    2.021075    2.704459
6.000000e+01    2.000298    2.660283
1.200000e+02    1.979930    2.617421
1.000000e+10    1.959964    2.575829

2.9 通過テスト

与えられた情報をそのままパラメータに渡せるので、Rのときより「ちゃんとやってる」感がある。別にだからどうしたということは無いのだが。

t.interval(alpha=0.95, df=499, loc=65, scale=np.sqrt(60)/np.sqrt(500))