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

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