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))

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

ハンバーガー統計学にようこそ!の第1章をPythonでやってみたもの。

準備

# -*- coding:utf-8 -*-
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

1.1 ポテトの長さの平均は

# ワクワクバーガー
waku = pd.Series(
         [3.5, 4.2, 4.9, 4.6, 2.8, 5.6, 4.2, 4.9, 4.4, 3.7,
          3.8, 4.0, 5.2, 3.9, 5.6, 5.3, 5.0, 4.7, 4.0, 3.1,
          5.8, 3.6, 6.0, 4.2, 5.7, 3.9, 4.7, 5.3, 5.5, 4.7,
          6.4, 3.8, 3.9, 4.2, 5.1, 5.1, 4.1, 3.6, 4.2, 5.0,
          4.2, 5.2, 5.3, 6.4, 4.4, 3.6, 3.7, 4.2, 4.8]
          )

# モグモグバーガー
mogu = pd.Series(
        [4.5, 4.2, 3.9, 6.6, 0.8, 5.6, 3.2, 6.9, 4.4, 4.7,
        3.8, 3.0, 3.2, 4.9, 7.6, 3.3, 7.0, 3.7, 3.0, 4.1,
        5.8, 4.6, 4.0, 2.2, 7.7, 3.9, 6.7, 3.3, 7.5, 2.7,
        5.4, 5.8, 5.9, 3.2, 5.1, 3.1, 6.1, 4.6, 2.2, 4.0,
        6.4, 5.2, 3.3, 6.4, 6.4, 2.6, 2.6, 5.2, 5.8]
        )

# 集計
# 総和と平均を求める
waku.sum() # 224.0
waku.mean() # 4.5714285714285712
mogu.sum() # 226.10000000000002
mogu.mean() # 4.6142857142857148

箱ひげ図はDataFrameにしてから描画すると簡単。

df = pd.DataFrame({"waku":waku, "mogu":mogu})
df.boxplot(grid=False)
plt.show()

1.2 度数分布

度数分布はpandas.cutで分割し、これに基づいてpandas.groupbyで元データを分割、グループごとにcountを適用するという手順をとることで得られる(cf. pandas で年齢階級をつくる - Qiita)。

waku.groupby(pd.cut(waku, np.arange(0, 9, 1))).count()
(0, 1]     0
(1, 2]     0
(2, 3]     1
(3, 4]    14
(4, 5]    19
(5, 6]    13
(6, 7]     2
(7, 8]     0
dtype: int64
mogu.groupby(pd.cut(mogu, np.arange(0, 9, 1))).count()
(0, 1]     1
(1, 2]     0
(2, 3]     7
(3, 4]    13
(4, 5]     8
(5, 6]     9
(6, 7]     8
(7, 8]     3
dtype: int64

ヒストグラムの方が作るのは簡単な印象。

df.hist(bins=np.arange(0,9,1),
        sharey=True,  # y軸を揃える
        grid=False,   # グリッド非表示
        layout=(2,1)) # 2行1列に並べる
plt.show()

1.4 分散と標準偏差

df.var(ddof=0) # 分散
mogu    2.584898
waku    0.680816
dtype: float64
df.std(ddof=0) # 標準偏差
mogu    1.607762
waku    0.825116
dtype: float64

ddofDelta Degrees of Freedomで、自由度はサンプルサイズをNとしてN - ddofとして指定される。Pandasの場合はデフォルトが1なので、何も指定しないと不偏分散およびそれに基づく標準偏差が計算される。今回はテキストに合わせるためにあえて0を指定している。

なお、numpy.varおよびnumpy.stdも分散および標準偏差を計算するメソッドであり、ddofという同名の引数が用意されているが、なんとこちらはデフォルトが0である。

1.9 通過テスト

df = pd.DataFrame({
  "桜組" : [78, 62, 81, 59, 72, 68, 75, 65, 80, 60, 78, 62, 70],
  "桃組" : [70, 72, 68, 75, 65, 71, 69, 76, 64, 80, 60, 73, 67],
  "柳組" : [57, 59, 55, 62, 52, 58, 56, 63, 51, 67, 47, 60, 54]
  })

平均、分散、標準偏差df.mean(), df.var(), df.std()でそれぞれ計算できるが、df.describe()を使うとデータ数、平均、標準偏差といわゆる五数要約(最大値、75%点、中央値、25%点、最小値)を確認できる。

df.describe()
柳組   桃組  桜組
count   13.000000   13.000000   13.000000
mean    57.000000   70.000000   70.000000
std 5.400617    5.400617    7.937254
min 47.000000   60.000000   59.000000
25% 54.000000   67.000000   62.000000
50% 57.000000   70.000000   70.000000
75% 60.000000   73.000000   78.000000
max 67.000000   80.000000   81.000000

また、scipy.stats.describe()ではデータ数、最大値、最小値、平均、分散、歪度、尖度を確認でき、ddofで自由度を指定できる。

import scipy as sp
sp.stats.describe(df)
DescribeResult(nobs=13, minmax=(array([47, 60, 59]), array([67, 80, 81])), mean=array([ 57.,  70.,  70.]), variance=array([ 29.16666667,  29.16666667,  63.        ]), skewness=array([ 0.,  0.,  0.]), kurtosis=array([-0.44902857, -0.44902857, -1.47721928]))

ついでに日本語ラベル付きのヒストグラムでも描いてみようと思って日本語表示に結構手間取った。seabornの設定でフォントを指定するとmatplotlibの設定もまとめて上書きされるので楽っぽい(cf.matplotlib と Seaborn の軸の日本語設定 - Qiita)。

import seaborn as sns
sns.set(font='Ricty') # Rictyは別途インストールしたフォント

df.hist()
plt.show()

seabornでも描いてみよう。

for i in df.columns:
  ax = sns.distplot(df[i],
                    label=i,
                    bins = 4)

ax.legend()
plt.show()

適当にやっても綺麗に出力される印象(cf.python - Two seaborn distplots one same axis - Stack Overflow)。

Rによるハンバーガー統計学 #3

(引き続きハンバーガー統計学にようこそ!の第3章をRでやってみたものです。)

3.1 チキンの売上は少ないのか

このようなタイプのデータを「とりあえずプロットする」方法として、モザイクプロットというものがある。

sales <- matrix(c(435, 265, 165, 135), nrow = 2)
rownames(sales) <- c("wakuwaku", "mogumogu")
colnames(sales) <- c("potato", "chicken")
mosaicplot(sales, col=TRUE) 

f:id:Rion778:20170802220451p:plain

モグモグのチキン売上割合が高そうなのがなんとなく見て取れるだろう。

もし差が無かったとした場合の度数(期待度数)は、ちょっとズルいやり方だけど、chisq.test()の中の$expectedを覗けば分かる。

> (expected <- chisq.test(sales)$expected)
         potato chicken
wakuwaku    420     180
mogumogu    280     120

3.2 カイ二乗値とカイ二乗分布

カイ二乗値を計算してみよう。

> dif <- sales - expected # 期待度数との差を求める
> mean(dif)  # 平均はゼロ😌
[1] 0
> sum(dif^2) # 自乗してから足せば…🤔
[1] 900
> sum(dif^2 / expected) # 期待度数で割ってから足せば…🤔🤔
[1] 4.464286

ピンポン玉実験の部分もやってみよう。

> x <- 5:9
> (x - 5)^2 / 5 + (10 - x - 5)^2 / 5
[1] 0.0 0.4 1.6 3.6 6.4

自由度1のカイ二乗分布をプロットしてみよう。

chisq <- seq(0, 8, by = 0.01)
plot(chisq, dchisq(chisq, df = 1), 
     type = "l",
     ylim = c(0, 1),
     xlab = expression(chi^2),
     ylab = "p")

f:id:Rion778:20170802221255p:plain

なおこの例、カイ二乗分布から計算されるp値は近似値なので注意。「オレンジ50:白50の箱から10個ピンポン玉を取り出す」という事象は、箱の外と中のピンポン玉の数を調べることで2x2の分割表で表現でき、従って対応するp値は厳密に計算できる。しかもセルの中には期待度数が5というものが含まれるので、結構ずれるんじゃないかと思われる。

先程描画した自由度1のカイ二乗分布の上に、先程の定義で計算されたカイ二乗値とフィッシャーの正確確率検定で計算されるp値の対応をポイントで示してみよう。

fisher_p <- function(N, ...){
  sapply(N, 
         function(n){
            fisher.test(matrix(c(50-n, n, 50-(10-n), 10-n), nrow=2), ...)$p.value
         }
  )
}
chisq_point <- (x - 5)^2 / 5 + (10 - x - 5)^2 / 5
points(chisq_point, fisher_p(5:9))

f:id:Rion778:20170802231028p:plain

では次に自由度別のカイ二乗分布を示そう(グラフの描き方は正規分布など好きな関数の曲線を描きたい - Qiitaが参考になった)。

df = c(1, 3, 5)
ggplot(data.frame(x=c(0, 8)), aes(x=x)) +
  mapply(function(d, co){stat_function(fun=dchisq, args=list(df=d), aes_q(color=co))},
         df, factor(df)) +
  scale_color_discrete(name = "df") +
  labs(x = expression(chi^2),
       y = "p") +
  theme_classic() +
  theme(legend.position = c(0.8, 0.9),
        legend.direction = "horizontal")

f:id:Rion778:20170802231300p:plain

自由度が大きいというのは、分割表で言えばセルの数が多いということに相当する。セルの数が増えれば、それだけ期待値からの相対的なずれの合計も大きくなりやすいし、ずれのばらつきも大きくなるということが見て取れる。

3.3 カイ二乗検定

最初に計算したカイ二乗値と、自由度1のカイ二乗分布の1%点と5%点を比較してみよう。

> sum(dif^2/expected) # そもそもカイ二乗値はいくつだっけ…
[1] 4.464286
> qchisq(p=0.95, df = 1) # 左から積分していって0.95になる点
[1] 3.841459
> qchisq(p=0.99, df = 1)
[1] 6.634897

5%水準なら帰無仮説を棄却できるし、1%ならできないということが分かった。

が、もちろん最初の分割表をそのままchisq.test()に放り込んでいい。

> chisq.test(sales)

    Pearson's Chi-squared test with Yates' continuity correction

data:  sales
X-squared = 4.1716, df = 1, p-value = 0.04111

今回は期待度数が大きいのでカイ二乗分布で近似してもさほど支障はないが、fisher.test()によってフィッシャーの正確確率検定を実施しても良い。

> fisher.test(sales)

    Fisher's Exact Test for Count Data

data:  sales
p-value = 0.041
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.010939 1.782292
sample estimates:
odds ratio 
  1.342623 

ハンバーガーを入れてみる

唐突なタイトル回収だ。

> sales2 <- cbind(sales, hamburger = c(650, 350)) # 唐突なタイトル回収
> sales2
         potato chicken hamburger
wakuwaku    435     165       650
mogumogu    265     135       350
> mosaicplot(sales2, col = TRUE)

f:id:Rion778:20170802232123p:plain

> chisq.test(sales2)

    Pearson's Chi-squared test

data:  sales2
X-squared = 9.9048, df = 2, p-value = 0.007067

> fisher.test(sales2)

    Fisher's Exact Test for Count Data

data:  sales2
p-value = 0.007509
alternative hypothesis: two.sided

3.9 通過テスト

> d <- matrix(c(24, 8, 18, 18), nrow = 2)
> rownames(d) <- c("kokugo", "sansu")
> colnames(d) <- c("sakura", "momo")
> d
       sakura momo
kokugo     24   18
sansu       8   18
> mosaicplot(d, col=TRUE)

f:id:Rion778:20170802232937p:plain

この問題は有意水準1%だと帰無仮説は棄却できず、有意水準5%だとカイ二乗検定かフィッシャーの正確確率検定かで判断が異なる例となっている。そんな些細な差で判断を分けるよりおとなしく追試をしたほうが良いだろうが。

> chisq.test(d)

    Pearson's Chi-squared test with Yates' continuity correction

data:  d
X-squared = 3.4874, df = 1, p-value = 0.06184

> fisher.test(d)

    Fisher's Exact Test for Count Data

data:  d
p-value = 0.04643
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.9600903 9.7553314
sample estimates:
odds ratio 
  2.950439 

モザイクプロットでは結構差があるように見えるのにp値がそこそこ大きくて、カイ二乗検定とフィッシャーの正確確率検定で結果が異なるのは、要するにサンプルサイズが小さいのだ。

それよりどうもこの問題、ここまでのテキストを参照するに「帰無仮説を採択する」と答えさせようとしているように思えるのだが、帰無仮説は採択できないのでその点には注意してもらいたい(cf. 「帰無仮説を採択」? | Okumura’s Blog)。

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

(引き続きハンバーガー統計学にようこそ!の第2章をRでやってみたものです。)

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

前回も説明したとおり、標準のvar()sd()でそれぞれ不偏分散、標本標準偏差が計算される。

wakupote <- c(47, 51, 49, 50, 49, 46, 51, 48, 52, 49)
mean(wakupote); var(wakupote); sd(wakupote)

2.3 信頼区間/区間推定

テキスト通りにやるならこのように。

# t値を求める
qt(0.025, 9) # 自由度9の場合の95%信頼区間の下側t値 qはquantileから
qt(0.975, 9) # 同上側t値
qt(0.005, 9) # 自由度9の場合の99%信頼区間の下側t値
qt(0.995, 9) # 同上側t値

# 信頼区間を求める
se <- function(x){ sd(x) / sqrt(length(x)) } # 標本標準誤差を求める関数

# 95%信頼区間
mean(wakupote) + qt(0.025, 9) * se(wakupote) # 下側
mean(wakupote) + qt(0.975, 9) * se(wakupote) # 上側
# 99%信頼区間
mean(wakupote) + qt(0.005, 9) * se(wakupote) # 下側
mean(wakupote) + qt(0.995, 9) * se(wakupote) # 上側

もちろん、Rを使っているのであればt.test()に突っ込むだけという姿勢が正しい。必要であれば信頼区間だけを取り出すこともできる。

> t.test(wakupote, conf.level = 0.95)$conf.int
[1] 47.85957 50.54043
attr(,"conf.level")
[1] 0.95
> t.test(wakupote, conf.level = 0.99)$conf.int
[1] 47.27432 51.12568
attr(,"conf.level")
[1] 0.99

t値が計算できるのだから実用的な意味は特に無いがt分布表も作ってみよう。

t_table <- data.frame(
  df = (df <- c(1:30, 40, 60, 120, Inf)),
  p95 = qt(0.975, df),
  p99 = qt(0.995, df)
)
knitr::kable(t_table)
df p95 p99
1 12.706205 63.656741
2 4.302653 9.924843
3 3.182446 5.840909
4 2.776445 4.604095
5 2.570582 4.032143
6 2.446912 3.707428
7 2.364624 3.499483
8 2.306004 3.355387
9 2.262157 3.249836
10 2.228139 3.169273
11 2.200985 3.105806
12 2.178813 3.054540
13 2.160369 3.012276
14 2.144787 2.976843
15 2.131449 2.946713
16 2.119905 2.920782
17 2.109816 2.898230
18 2.100922 2.878440
19 2.093024 2.860935
20 2.085963 2.845340
21 2.079614 2.831360
22 2.073873 2.818756
23 2.068658 2.807336
24 2.063899 2.796940
25 2.059539 2.787436
26 2.055529 2.778714
27 2.051830 2.770683
28 2.048407 2.763263
29 2.045230 2.756386
30 2.042273 2.749996
40 2.021075 2.704459
60 2.000298 2.660283
120 1.979930 2.617421
Inf 1.959964 2.575829

2.4 平均的チキンを推定する

🐔

chicken <- c(568, 530, 581, 554, 536, 518, 564, 552)
mean(chicken)
var(chicken)
sd(chicken)
se(chicken)
t.test(chicken, conf.level = 0.95)$conf.int
t.test(chicken, conf.level = 0.99)$conf.int

2.9 通過テスト

  • (1)
    1. 母集団
    2. 標本
    3. 無作為抽出
    4. 不明である

dについて実際どのような解答が想定されているか分からないが、問題文の情報量だけで度数分布の形状を推定することは不可能である。例えば、テストの総問題数が1問であった場合を考えてみれば想像できるだろう。実際にはどのような分布も想定すべきでないという姿勢で分析にかかるべきではないのか。

  • (2)

もっとうまいやり方があるような気がして仕方がない。

# 95%信頼区間
65 + qt(0.025, 500) * sqrt(60)/sqrt(500)
65 + qt(0.975, 500) * sqrt(60)/sqrt(500)
# 99%信頼区間
65 + qt(0.005, 500) * sqrt(60)/sqrt(500)
65 + qt(0.995, 500) * sqrt(60)/sqrt(500)
  • (3)

信頼区間の計算式においては確率pのみが変数であるので、逆に計算された信頼区間を変数として確率pが計算されるような関数を考えることができる。このとき、pの範囲が2.5%〜97.5%となるような区間が95%信頼区間、同様に0.5%〜99.5%となるような区間が99%信頼区間である。

「信頼区間に母平均が含まれる確率が95%または99%である」という説明は信頼区間の説明としては正しくない(例えばRで楽しむ統計 (Wonderful R 1)に詳しい解説がある。上述の方法とは異なる信頼区間の構成方法も紹介されている。)。

Rによるハンバーガー統計学 #1

(ハンバーガー統計学にようこそ!をRでやってみたものです。)

1.1 ポテトの長さの平均は

# ワクワクバーガーのポテトの長さデータ
waku <- c(3.5, 4.2, 4.9, 4.6, 2.8, 5.6, 4.2, 4.9, 4.4, 3.7,
          3.8, 4.0, 5.2, 3.9, 5.6, 5.3, 5.0, 4.7, 4.0, 3.1,
          5.8, 3.6, 6.0, 4.2, 5.7, 3.9, 4.7, 5.3, 5.5, 4.7,
          6.4, 3.8, 3.9, 4.2, 5.1, 5.1, 4.1, 3.6, 4.2, 5.0,
          4.2, 5.2, 5.3, 6.4, 4.4, 3.6, 3.7, 4.2, 4.8)

# 総和を求める
sum(waku)

# 平均を求める
sum(waku)/49

# 平均を関数で求める
mean(waku)

# モグモグバーガーのポテトの長さデータ
mogu <- c(4.5, 4.2, 3.9, 6.6, 0.8, 5.6, 3.2, 6.9, 4.4, 4.7,
          3.8, 3.0, 3.2, 4.9, 7.6, 3.3, 7.0, 3.7, 3.0, 4.1,
          5.8, 4.6, 4.0, 2.2, 7.7, 3.9, 6.7, 3.3, 7.5, 2.7,
          5.4, 5.8, 5.9, 3.2, 5.1, 3.1, 6.1, 4.6, 2.2, 4.0, 
          6.4, 5.2, 3.3, 6.4, 6.4, 2.6, 2.6, 5.2, 5.8)

# 平均値を求める
mean(mogu)

1.2 度数分布

実はhist()の戻り値の中に度数分布表は含まれていたりするのだけれど、このタイミングでそれもどうかなと思うので別解。

table(cut(waku, breaks = 0:8))
table(cut(mogu, breaks = 0:8))

cut()numericを指定の区間で区切ったfactorに変換する。これをtable()で集計すれば度数分布表が得られる。結果を見たほうが早いだろう。

> table(cut(waku, breaks = 0:8))

(0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] 
    0     0     1    14    19    13     2     0 
> table(cut(mogu, breaks = 0:8))

(0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] 
    1     0     7    13     8     9     8     3 

1.4 分散と標準偏差

標準のvar()sd()は不偏分散に基づいた標本分散と標本標準偏差を計算するので、テキストに沿った値が欲しければ関数を定義する必要がある。

my_var <- function(x){mean((x - mean(x))^2)}
my_sd  <- function(x){sqrt(my_var(x))}
# 組み込みのvar(), sd()は不偏分散に基づくので計算結果が異なる

# ワクワクバーガーの分散と標準偏差
my_var(waku)
my_sd(waku)

# モグモグバーガーの集計
my_var(mogu)
my_sd(mogu)

1.5 チキンで行こう

waku_c <- c(135, 142, 149, 146, 149, 144, 136, 138,
            156, 153, 150, 147, 136, 160, 142, 157)
mogu_c <- c(144, 143, 139, 166, 169, 144, 147, 138, 
            176, 133, 170, 137, 146, 140, 122, 177)
# 平均、分散、標準偏差の計算
mean(waku_c); my_var(waku_c); my_sd(waku_c)
mean(mogu_c); my_var(mogu_c); my_sd(mogu_c)

1.9 通過テスト

score <- data.frame(
  桜組 = c(78, 62, 81, 59, 72, 68, 75, 65, 80, 60, 78, 62, 70),
  桃組 = c(70, 72, 68, 75, 65, 71, 69, 76, 64, 80, 60, 73, 67),
  柳組 = c(57, 59, 55, 62, 52, 58, 56, 63, 51, 67, 47, 60, 54)
)

my_summary <- function(x){
  data.frame(
    平均     = round(mean(x), 2),
    分散     = round(my_var(x), 3),
    標準偏差 = round(my_sd(x), 3)
  )
}

apply(score, 2,  my_summary)

なんでn-1なんですか

不偏分散はなぜnではなくn-1で割るのか、という問題。

以前似たようなことやったけどきちんと証明していなかったのと今ならどう書くか試してみたかったので。

証明をする

確率変数Xの実現値をX_i、平均を\bar{X}とする。

まず、\sum\left(X_i - \bar{X}\right)^2の期待値を求めよう。

任意の数\muを足して引く。


E\left[\sum\left((X_i - \mu) - (\bar{X} - \mu)\right)^2\right]

ここで、


\sum(x_i - \bar{x})^2 = \sum x_i^2 - n\bar{x}^2

であることに注意すれば、


\sum E \left[ (X_i - \mu)^2 \right] - nE \left[ (\bar{X} - \mu)^2 \right]


\sum E \left[ (X_i - \mu)^2 \right] - nE \left[ \left( \frac{1}{n} \sum (X_i - \mu) \right)^2 \right]

ここで任意の数\muが母平均に等しかったとすれば、


\sum E \left[ (X_i - \mu)^2 \right] = \sigma^2

であるので、


n\sigma^2 - \sigma^2 = (n-1)\sigma^2

従って、\sigma^2の不偏推定量\hat{\sigma^2}


\hat{\sigma^2} = \frac{1}{n-1}\sum\left(X_i - \bar{X}\right)^2

となる。

証明をしない

標準正規分布からのサンプリングに対し、nで割った場合とn-1で割った場合の分散の分布をシミュレーションしてプロットしてみよう。

library(dplyr)
library(ggplot2)

# (n-1)/nを乗ずることでnで割った場合の分散を求める
varp <- function(x){
  var(x) * (length(x) - 1)/length(x)
}

# n = 2:20において各1000回ずつvar(rnorm(n))及びvarp(rnorm(n))を求める
d <- data.frame()
N <- 1000
for(n in 2:20){
  d <- rbind(d,
             data.frame(n = n,
                        type = "var",
                        var = replicate(N, var(rnorm(n)))))
  d <- rbind(d,
             data.frame(n = n+0.3, # プロット時にちょっと右にずらすため調整
                        type = "varp",
                        var = replicate(N, varp(rnorm(n)))))
}

ggplot(d, aes(x = n, y = var, color = type)) + 
  geom_point(alpha = 0.1) +
  stat_summary(fun.y = mean, geom = "line") + # 平均をプロット
  theme_classic() +
  geom_hline(yintercept = 1, lty = 2, col = "gray", alpha = 0.3) # 真の分散=1を示す線

f:id:Rion778:20170706233312p:plain

このようにnで割った分散の期待値は真の分散(=1)より少し小さくなる傾向がある。一方でn-1で割った場合の期待値は真の分散に一致する。また、nで割った場合でもnが大きくなれば真の分散に近づく傾向があることが分かる。実際、nで割った場合でもn→∞で真の分散に一致するので、nで割った分散は不偏推定量ではないが一致推定量ではある。

ユーザーの平均継続期間が解約率の逆数であることを等比数列を使わずに理解する

ユーザーの平均継続期間は「解約率=一定」と仮定すれば「1/解約率」で示せる。

これについてはユーザの平均継続期間が「1/解約率」で求められることの数学的証明 - it's an endless world.等比数列を用いたわかりやすい証明があるが、「解約率=一定」とした場合の平均継続期間は「故障率=一定」とした場合のMTBFと計算上は同じなので、MTBFの計算に関する考え方を流用すると等比数列を使わずに示すこともできる。

MTBF(平均故障間隔)

MTBF(Mean Time Between Failures = 平均故障間隔)はあるシステムが故障するまでの平均動作時間で、システムの稼働時間/故障回数として求める。

例えば100個の装置を100時間観察し、その間に2つ装置が故障した場合、MTBF


\mathrm{MTBF} = \frac{100 \times 100}{2} = 5,000\ (\mathrm{hours})

のように延べ稼働時間から求める(cf. 平均故障間隔 - Wikipedia)。要するに、100個の装置100時間稼働を、1個の装置10000時間稼働と同一視する。

ユーザーの平均継続期間

現在のユーザー数をx_tdt時間経過後のユーザー数を x_{t+dt}とすれば、その間の解約ユーザー数はx_t - x_{t+dt} = -dx_tである*1。このとき、平均継続期間Uは、次式で求められる。

U = \frac{x_t dt}{-dx_t}

要するに、ユーザー1人についてx_t dtという期間中に-dx_t回解約があるとみなして、平均契約期間を計算する。サービスを解約してもなぜかすぐに再開する1人のユーザーを仮定するのだ。これは上述のMTBF計算の場合と同様の考え方である。

次に、単位時間中のユーザーの離脱率が一定値rだとすると、ユーザー数は現在のユーザー数x_tに離脱率rを乗じた速度で減少する。

すなわち、ユーザー数の変化速度\frac{dx_t}{dt}は次式で表される。

\frac{dx_t}{dt} = -rx_t

これを先程の式と比較すれば、

U = \frac{1}{r}

であり、平均継続期間は離脱率の逆数であることが示せる。

*1:dt = (t+dt) - tdx_t = x_{t+dt} - x_tであるため、dx_tx_tが時間とともに減少するという仮定のもとでは負の値であり、これにマイナスを付けたものが解約ユーザー数となる

Project Euler Problem 22

Problem 22 - PukiWiki

  • 名前のリストが入ったテキストファイルを読み込み、ソートする。
  • アルファベットを1から26までの整数に対応付けてその和を計算し、さらに名前の順位を乗ずることで「名前のスコア」を算出する。
  • 例えばCOLINなら3 + 15 + 12 + 9 + 14 = 53で、順位は938なので、53*938=49714がCOLINのスコアになる。
  • 全ての名前のスコアの和を求めよ。

解答

データの読み込みと基本的な文字列操作と演算ができれば特に問題はない…。

続きを読む

Project Euler Problem 21

問題の概要

Problem 21 - PukiWiki

  • d(n)をnの「真の約数」の和とする。
  • 「真の約数」はnの約数のうちnを除外したものと定義する。
  • d(a) == b かつ d(a) == bであるようなaとbをAmicable numbers(友愛数)と呼ぶ。
  • 完全数(すなわちa == b)は除外する。
  • 10000以下の友愛数の総和を求めよ。

解答

最初完全数を除外する条件を見逃していてハマった。

ちなみに問題文からはaが10000以下でbが10001以上の場合の処理方法が分からないが、とりあえず今回の範囲ではそのような友愛数のペアは出現しない。

続きを読む

スペインの施設園芸

スペインにはハウスがたくさんある

昔作った資料をもとに話す機会があったので、スペインの施設園芸について手持ちの文献を再度調べていた。

スペインの南部のアルメリア県にはハウスが集中している地域がいくつかあって、特にアルメリア市の西部に位置するアドラ市とエルエヒド市を中心とするエリアは明確にハウスが密集している。衛星写真に切り替えて確認してみてほしい。

このエリアの施設面積だけで2万haあると言われており、日本全国の施設面積の1/3程度に相当する。

スペインの施設はブドウの棚から発展した

この地域でこのように施設園芸が普及するまでの経緯が面白い。もともとこの地域はブドウ(ワイン用ではなく生食用だったそうだ)の栽培が盛んで、棚仕立てでブドウを栽培していたのだが、ある時生産者がこの棚の上にビニールを載せることでその下で色々な野菜を作れる、ということに気付いたそうだ。

そこから施設園芸が広まったので、ハウスの形状も独特であった。今もこの地域で最も普及しているハウスはRaspa y Amagado(ラスパ・イ・アマガド)と呼ばれるもので、ベースはビニールの載ったブドウの棚なのだが、雨水を排水するために屋根がゆるくギザギザと波打っている。Raspaは魚の骨、Amagadoは屈曲を意味し、天井部の高い部分がRaspaで谷になっている部分がAmagadoなのだそうだ。

この話を「施設と園芸」の2012年の夏号(No.158)で読んで、一度見てみたいものだ(写真は載っていたがモノクロで分かりにくかった)と思っていたのだが、昨日改めて航空写真などを確認していてそう言えばストリートビューというものがあったと気付き、適当な場所を確認してみたら記事の説明そのもののRaspa y Amagadoを確認することができた。天井部の高い部分はまさに背骨という感じだ。

上に示した場所の場合は最新の写真が2011年7月のものだが、ハウスの中は空っぽだ。この地域では夏は生産物単価が下がるのと、ハウス内部が高温になって作業が大変という理由から基本的に9月から翌年5月までの作型で栽培をするそうだ。同じ場所の2009年1月の写真を遡って確認できるので見てみると、ハウスの中にトマトが植わっているのが分かる。

夏の間は前作の片付けや次作の準備を進めつつビーチで優雅に過ごすそうだ。

ハウスはなんだか安っぽく見える(実際建造費は安いのかもしれない)が、開口部はかなり目合の細かいネットで完全に覆われているのが伺える。出入り口は基本的に二重扉となっており、外部からの害虫の侵入を防止している。中には送風装置を備えて体に付着した害虫を払い落とすような仕組みを備えた施設もあるそうだ。日本の先進地域並の装備である。

加えてスペインの農業で特徴的なのは天敵の利用が一般的であるという点だ。天敵利用は高度な技術を必要とするが、ヨーロッパでは有機農業への需要が高いこともあり、天敵利用技術が発達しており、天敵などの生物農薬を生産販売するメーカーも複数ある。販売する苗にあらかじめ天敵を定着しておくということも行われているそうだが、これは生産者の天敵利用に対する深い理解と技術力がなければ到底不可能なことだ。天敵利用に関しては明らかに日本の先を行っている。

かつては気候に恵まれ安く野菜が生産できるという点だけが利点で、オランダのような施設園芸の先進国は天敵利用技術に先んじていたため、多少値段が高くても有意性を維持できていた。というようなことが古い教科書には書いてある。

そのスペインに天敵利用が急速に広まり、HACCPやグローバルGAPを取得する生産者も続々出てきて施設園芸の先進国も苦しい立場に…というのも少し昔の話で、今ではモロッコやイスラエルが農業新興国として安くて高品質な野菜を輸出するようになっており、今度はスペインの立場も危うくなっているそうだ。

などということが書いてあった記事に関心していたのももう5年も前の話だ。5年もたてば状況も変わっている可能性が高い。また時間をみつけて新しい情報を確認したいところだ。