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年もたてば状況も変わっている可能性が高い。また時間をみつけて新しい情報を確認したいところだ。