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

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

以下解答。

2017-08-06

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

Python Statistics

引き続きハンバーガー統計学にようこそ！の第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))

2017-08-06

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

Python Statistics

ハンバーガー統計学にようこそ！の第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

ddofはDelta 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)。

2017-08-02

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

R Statistics

(引き続きハンバーガー統計学にようこそ！の第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）。

2017-07-20

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

R Statistics

(引き続きハンバーガー統計学にようこそ！の第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問であった場合を考えてみれば想像できるだろう。実際にはどのような分布も想定すべきでないという姿勢で分析にかかるべきではないのか。

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

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

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

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

2017-07-20

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

R Statistics

(ハンバーガー統計学にようこそ！を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)

2017-07-06

なんでn-1なんですか

R Statistics

不偏分散はなぜ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で割った分散は不偏推定量ではないが一致推定量ではある。

2017-06-06

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

Math

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

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