Cohen の d

中澤さんが

> effsizeパッケージのcohen.d()関数の計算は何か変だ。大久保・岡田(2012)のSpの式で標本分散のはずのところが，不偏分散になっている。どちらかがCohenの定義を誤解していると思われるが，不偏分散にサンプルサイズを掛けるのは筋が通らないので，たぶんeffsizeパッケージが間違っているのだろう。

と書いている。

大久保・岡田(2012) がないのでよく分からないのだが，Cohen's d は，奥村先生のページにもあるが，

で，

stdev = sqrt(((n1 - 1) * s[1]^2 + (n2 - 1) * s[2]^2)/(n1 + n2 - 2))

だからいいのではないか？

effsize は

Version:            0.7.6
Date:               2019-07-17

それとも，私が何か勘違いしている？

追記：中澤さんの勘違いだったそうだ。この記事のあと，修正が入ったようだ。

> x = c(1,2,3,2,3,4)
> y = c(5,4,3,4,5,6,5)
> cohen.d(c(x, y), factor(rep(1:2, c(6, 7))))

Cohen's d

d estimate: -2.051542 (large)
95 percent confidence interval:
     lower      upper
-3.5627102 -0.5403745 

> mes(mean(x), mean(y), sd(x), sd(y), n1, n2)
Mean Differences ES:
 
 d [ 95 %CI] = -2.05 [ -3.56 , -0.54 ]
  var(d) = 0.47
  p-value(d) = 0.01
  U3(d) = 2.01 %
  CLES(d) = 7.34 %
  Cliff's Delta = -0.85
 
 g [ 95 %CI] = -1.91 [ -3.31 , -0.5 ]
  var(g) = 0.41
  p-value(g) = 0.01
  U3(g) = 2.82 %
  CLES(g) = 8.86 %

ぼったくり宝くじ

オイ，録音してないやろな

なんですか，いきなり。してませんがな。なんでせなならんの。

ちょっとヤバい話だ。良く効け。おれらも相当あくどい仕事してるけど，世の中には堂々とひどい商売やってる奴がいるぜ。

なんですか？

宝くじや。

あ，今，「買わないという選択はない」とかコマーシャルで言ってる。

そや。

でも兄貴。あのお堅い銀行さんも噛んでいる，まともな商売じゃないんですか。

どこがまともじゃ。アホンダラ。ええか，100万円分の宝くじ買ったとするやろ，なんぼ当たるかわかるか？

え〜っ。100 万円も買うんでっか。1000万位は当たりそう。

あほ抜かせ。ええとこ行って，46 万位しか当たらんのじゃ。

え〜っ。え〜っ。なんでそんな。

ええか。ちゃんと，総務省，総務省やぞ，の Web ページというか PDF

http://www.soumu.go.jp/main_content/000084191.pdf

に，ちゃんと書いてある。

実行還元率の所に，45.7% と書いてあるやろが。

わ。ビックリした。かたぎの商売とは思えんピンハネ，ぼったくり。

宝くじは，買うたんびに，半分以下になるんや。100万円からスタートしたら，何回目で宝くじの一枚も買えない金額 300 円になるかわかるか？

いや，わての頭にはコンピュータ入ってませんから，わかりませんわ。

ほなら，ちょっと計算したる。

> a = 1000000
> for (i in 1:11) {
+     a = round(a * 0.457)
+     cat(sprintf("%2i, %6.0f\n", i, a))
+ }
 1, 457000
 2, 208849
 3,  95444
 4,  43618
 5,  19933
 6,   9109
 7,   4163
 8,   1902
 9,    869
10,    397
11,    181

な。途中 300 円未満の端数は無視するとして，11 回目には 181 円や。

一つ前の話と前の話は，ゆっくり減ったり増えたりする話やったけど，これはスゴイ。地獄の底へ真っ逆さまや。

ええか，宝くじには手を出したらあかんで。

へい。でも，見せてくれた図の右にある公営競技ってのが，気になりますね。還元率が74.8で，実行還元率でも58.5% じゃないですか。宝くじ買うより競馬や競輪やれって言うことですか。

数字の罠に気をつけろ。公営ギャンブルは一日に何回かある。例えば 5 回あるとするやろ。そしたら 2 日いって，全部賭けるとしたら，

> 0.585^10*1000000
[1] 4694.168

ええか，100 万からスタートして2日目の終わりには 4 千円ちょっとしか残らんのやぞ。

もし 100 万あったら，宝くじ買ったりギャンブルするより，きれいなチャンネーのいるところで贅沢したほうがええ，ゆうことですね。

そや。昔っから，お上はあの手この手で庶民の懐から金をふんだくっていくもんだと決まってる。ギャンブルは一切するな。近頃じゃ IR とかなんとかいって，また金をふんだくる算段をしているみたいだぜ。

おれたちよりあくどいね。

おれたちっていっても，反社じゃないよ。

ちゃんちゃん。

貯蓄額倍増？

ご隠居さん，コンチハ。

ああ，はっつあんか。まあお上がり。

ちょっとご隠居さんにうかがいたいことがあるんですけど。

なんでも聞いとくれ。

「定年後に2000万円必要だ！」なんて言ってる人がいるんだけど，ホントですか。

まあ，ほんとのようなうそのような。なんだろね。

今，銀行の利率なんて雀も気の毒がるほど低いじゃないですか。

そうだねえ，昔は年利率 5% なんてときもあったんだが，今じゃ 0.2% なんてことも普通だねえ。

あっしが，今 100 万円持ってるとして，年利 0.2% の複利で預けるとして，元利合計が 200 万円になるのは何年後です？

おや，驚くねえ。100 万も持ってるのかい？

冗談言っちゃいけね〜。持ってるとしてといったでしょうに。

そうさなあ，年利 0.2% なら 350 年ぐらい掛かるな。

ご冗談でしょ。その前にテメーの命がなくならあ。

そうさなあ。

でも，ずいぶん早く答えが出ましたね。ご隠居さんの頭にはコンピュータでも入ってるんですかい。

いやいやそうではない。70 を年利のパーセント数で割ってやれば，二倍になるまでの年数がわかるんだよ。

年利を p パーセントとする，2倍になる年数を n として
(1 + p/100) ^ n = 2 を解く
n * log(1 + p/100) = log(2)
n = log(2) / log(1 + p/100)
  = log(2) / (p/100 - (p/100)^2/2 + (p/100)^3/3 - (p/100)^4/4 + ...)
  ≒ 0.6931472 / (p/100)
  ≒ 70 / p

70 割る 0.2 ですか。小数で割り算は，あっしの頭じゃ無理だ。

700 を 2 で割ってもいいんだよ。

はあ，それで 350 年か。

65 歳で定年のときに退職金を 1000 万円もらっらとして，年利 5% で預けても，2000 万円になるのは 14 年かかるな。79 際まで，1000 万円に手を付けなければということだがな。

所詮無理ですね。あっしなんか，退職金なんかないんですから。どうすりゃいいんです？

まあ，生活を切り詰めて，2000 万円など要らない生活をするんだね。

じゃあ，今と同じじゃないですか。

そうそう，高望みしないのが一番じゃ。ほっほっほ。

笑い事じゃねえや。

 

土用ウナギの日

土用の丑の日に限らないが，よく食の老舗で「江戸時代から継ぎ足し継ぎ足ししてきたタレ」とかいう話があり，そのたびに，「ばかか」と思う

まあ，ウナギのタレがどういう単位で存在するのかはさておき，たとえば単位不明で n = 1000000 としよう。毎日毎日そのうちの千分の1である m = 1000 を使い，そのぶんだけ新しいものを補充するとしよう。

当初，n 個の要素は全て 1 を持っている。a = rep(1, n)

毎日毎日そのうちの r = sample(n, m, replace=TRUE) の要素が新しいものに置き換えられる。

つまり，使われた要素は 1 から 0 に変えられる。a[r] = 0

一日が終わったとき，当初の要素が幾つあるかは sum(a) でわかる。

というようなシミュレーションプログラムが以下のようになるであろう。

n = 1000000
m = 1000
a = rep(1, n)
day = 9999
conc = integer(day)
for (i in seq_len(day)) {
    r = sample(n, m, replace=TRUE)
    a[r] = 0
    t = sum(a)
    conc[i] = t
    cat(sprintf("day %3i, parts = %i\n", i, t))
}
plot(1:day, conc, type="l", xlab="day")

day   1, parts = 999000
day   2, parts = 998000
day   3, parts = 997003
day   4, parts = 996006
day   5, parts = 995011
　　：
day 1101, parts = 332353
day 1102, parts = 332005
day 1103, parts = 331680
day 1104, parts = 331357
day 1105, parts = 331023
　　：
day 9995, parts = 50
day 9996, parts = 50
day 9997, parts = 50
day 9998, parts = 50
day 9999, parts = 50

図にすると，以下のよう。見事に指数関数的減衰。

27年ほど経つと濃度は 0.00005 になる。

もし，一日に 1/100 を入れ替えるとすると，1000日（約2年9ヶ月）で 0.00005 になる。

江戸時代からだと，ほとんど1個の分子すら残っていまい。

まあ，「一晩の内に，古いタレが，新しいタレを教育して昔ながらの味になる」なんて，おとぎ話があるのかもしれないけど。

トリボナッチ数列

「再帰関数と非再帰関数」，「事象が連続発生するまでの試行回数」にも出てきたトリボナッチ数列だが，

「n 階段があるとき，1段昇る，2段昇る，3段昇るのどれかを選択して昇るとき何通りの昇り方があるか」という設問がされることもある。

よく，「アルゴリズムをわかりやすく表現できる再帰プログラム」の例としても出てくるのが以下のプログラム。フィボナッチ数列の拡張ということだ。

def step(n):
    if n == 0:
        return 1
    if n < 0:
        return 0
    return step(n-1) + step(n-2) + step(n-3)

でも，これは相当スタックも喰うし，実行時間も馬鹿にならない。

import time
s = time.time()
for i in range(30):
    print(i, step(i))
print(time.time()-s)

30 項計算するのに 46秒も掛かる。

これを再帰ではないプログラムにすると，（これもフィボナッチ数列の非再帰的求め方のプログラムだが）

s = time.time()
a, b, c, = 1, 1, 2
for i in range(3, 30):
    a , b, c = b, c, a + b + c
    print(i, c)
print(time.time()-s)

当たり前だが，わずか 0.0018 秒で終わる。

 

事象が連続発生するまでの試行回数

やあ，ジュン。元気かい？

元気だよ。

前に，ポケモン Go のことで教えてもらったけど，またいいかな？

いいよ。でも，ポケモン Go に興味がない人にも分かるような例で説明してくれる？

そうか〜。うんとね。あ，そうそう。野球にたとえるとこんなになるかな。「3 割バッターが，3 連続ヒットを打てるまでに何回打席に立たなければならないか」でいいかな。前に教えてもらった負の二項分布に似ている気もする。

ヒットを打ったら 1，打たなかったら 0 で表すと，打席の記録は 100110100111 のように最後の3桁が 111 になったら 3 連続ヒットを打ったということ。当然だけど，最後の 3 桁より左には 1 が 3 つ以上連続することはない。この例だと，12 打席目で 3 連続ヒットを達成ということだけど，確率変数 x は 3 連続ヒットまでの打席数として 9 と考えるほうがいいんだろうね。

そうだね。3 打席連続でヒットを打ったときは 111 で終わりなので，連続ヒットまでの打席数は 0 だから，そのほうが都合がいいよね。
打率を p としよう。x = 0 のときの確率はわかるよね。

それは，簡単。p * p* p = p^3 だね。p = 0.3 なら 0.0270 だ。

じゃあ，x = 1 のときは？

打席の記録は 0111 だね。初打席が 1 だと 1111 になるけど，4打席目はないので，x = 0 のときの話になる。
確率は，打席の記録が 0 のとき 0.7, 1 のとき 0.3 を順に掛けていけばよいから
0.7 * 0.3 * 0.3 * 0.3 = 0.7 * 0.3^3 = 0.0189 だ。

次は，x = 2 だけど，最後の5打席目を打つまでに 111 が出ないように注意しないといけないね。

打席の記録は 10111 か 00111 だから，それぞれの起きる確率は 0.7 * 0.3^4 = 0.00567，0.7^2 * 0.3^3 = 0.01323 で，併せて 0.0189。
x = 3 のときは，000111, 010111, 100111, 110111 の 4 通り。でも，0 と 1 の個数で整理して (0 の個数，1 の個数) のように表すと(3, 3) が 1 通り, (2, 4) が 2 通り, (1, 5) が 1 通りなので，確率は
 (0.7^3 * 0.3^3) + 2*(0.7^2 * 0.3^4) + (0.7 * 0.3^5) = 0.0189 になるね。
x = 1, 2, 3 のどの場合も確率が 0.189 というのもおもしろいね。

これを x = 4, 5, ... と間違いなくやっていくのも大変なので，プログラムを書こう。途中に 111 がある場合を除いて勘定するために打席記録の最後が 0111 で終わるものだけを対象にする。

library(e1071)
f = function(k, p = 0.3) {
  Z = NULL
  b = bincombinations(k)
  P = apply(b, 1, function(x) {
    if (! grepl("111", paste(x, collapse=""))) {
      ones = sum(x)
      zeros = k-ones
      Z <<- c(Z, ones)
      p^(ones+3) * (1-p)^(zeros+1)
    } else {
      0
    }
    })
  print(unname(table(Z)))
  return(sum(P))
}

という関数を作って，以下のように実行する。

> p2 = numeric(11)
> for (i in 1:11) {
+   p2[i] = f(i)
+ }
[1]   1   1
[1]   1   2   1
[1]   1   3   3
[1]   1   4   6   2
[1]   1   5  10   7   1
[1]   1   6  15  16   6
[1]   1   7  21  30  19   3
[1]   1   8  28  50  45  16   1
[1]   1   9  36  77  90  51  10
[1]   1  10  45 112 161 126  45   4
[1]   1  11  55 156 266 266 141  30   1

これを図1のようにまとめる。

図の赤字で示したところが今求めた部分だ。「"111" より前の桁数」x が 2 〜 12 の場合について求めている。

それぞれの場合の確率も求められているよ。

> cbind(p2)
              p2
 [1,] 0.01890000
 [2,] 0.01890000
 [3,] 0.01838970
 [4,] 0.01803249
 [5,] 0.01767528
 [6,] 0.01731807
 [7,] 0.01697050
 [8,] 0.01662969
 [9,] 0.01629563
[10,] 0.01596832
[11,] 0.01564757

x = 3 の行を見れば，上に書いたのと同じになっているね。
のとき，1, 2, 1 とあるのは，1 が ない場合が 1 通り，1 個ある場合が 2 通り，2 個ある場合が 1 通りということである。最後の 0111 も併せて 0, 1 は (3, 3), (2, 4), (1,5) なので 確率は (0.7^3 * 0.3^3) + 2*(0.7^2 * 0.3^4) + (0.7 * 0.3^5) = 0.0189 となっている。

ちなみに，図1の合計欄の 1,2,4,7,13,24,44,81,... を「オンライン整数列大辞典」で検索してみると「Tribonacci numbers」というものであることが分かる。1+2+4 = 7, 2+4+7 = 13, ... だね。この表自体がある規則性を持っているという傍証だね。

x がもっと大きい範囲だとどうなるんだろうか。

x が大きくなると前述のプログラムでは計算できないので，図1の「1〜x 桁にある 1 の個数」が 4 の列の数列，1, 6, 19, 45, 90, 161, ... を「オンライン整数列大辞典」で検索してみると，「(1+x+x^2)^n を展開したときの x^4 の係数」というのがぴったり。5, 6 の列もみてみようか。
それぞれの列は x^0, x^1, x^2, ... の係数の数列になっているみたいだね。

n が大きくなると (1+x+x^2)^n を展開するのは大変じゃない？

真っ正直に展開する必要はないよ。係数だけを書き出すと

n = 1 のときは

  1 1 1

n = 2 のときには

  1 1  1
    1  1  1
+      1  1  1
===============
  1  2  3  2  1

n = 3 のときには

  1  2  3  2  1
     1  2  3  2  1
+       1  2  3  2  1
======================
  1  3  6  7  6  3  1

とみたいに，単純に 1 桁ずつずらして 3 回たしざんするだけ。Excel で作ってみよう。

というわけで，図 2 ができたね。

でも，縦方向に見たときの数列は確かに図 1 と同じだけど，横方向に見ても図 1 には一致しないね。

よく回りを見回せば，図 1 の x = 11 の行は，図2の n=11，x^0 のマス目から右上へ黄色で色づけしたマス目の数値を読んでいけば得られる。すなわち，1, 10, 45, 112, 161, 126, 45, 4 だよ。

確率はだらだらと下がっていくようだけど，どんな風になるのかな。

この計算があっているかどうかも併せて，シミュレーションで確かめてみよう。ちょっと冗長なプログラムだけど，素直に書くと，

triple = function(p, k = 3) {
  count = -3
  while (TRUE) {
    if ({count = count + 1; runif(1) < p}) {
      if ({count = count + 1; runif(1) < p}) {
        if ({count = count + 1; runif(1) < p}) {
          return(count)
        }
      }
    }
  }
}

p = 0.3
x = replicate(1000000, triple(p))
plot(table(x))
median(x)
mean(x)

平均値は 48.54934，中央値は 33 ぐらいか。ずいぶん右裾が長いね。

計算結果ともよく合っていると思う。線が計算値，赤丸がシミュレーション値だよ。

x=12 のときと x=13 のときの確率が同じというのがちょっと気になるけど。

x  Probability Simulated
 0  0.02700000  0.026932
 1  0.01890000  0.018876
 2  0.01890000  0.018859
 3  0.01890000  0.018700
 4  0.01838970  0.018510
 5  0.01803249  0.018055
 6  0.01767528  0.017658
 7  0.01731807  0.017579
 8  0.01697050  0.017031
 9  0.01662294  0.016504
10  0.01629563  0.016123
11  0.01596832  0.015904
12  0.01564757  0.015522
13  0.01564757  0.015265
14  0.01533327  0.014866
15  0.01502529  0.014577
16  0.01472348  0.014570
17  0.01442774  0.014255
18  0.01413795  0.013888
19  0.01385397  0.013625
20  0.01357569  0.013273
21  0.01330301  0.013093

これですっきりした。ありがとう。

デカルトの正葉線

デカルトの正葉線 folium of Descartes を描いてみようと思い，やってみたが自力ではちょっと難しかった。

媒介変数表示だと，大抵は
x = 3*a*t/(1+t^3)
y = 3*a*t^2/(1+t^3)
が紹介されている。しかしこれだと，原点付近が描画できない。
そこで，更に
t = (1+s) / (1-s) と変数変換して，
x = 1.5*a*(1-s)*(1-s^2) / (1+3*s^2)
y = 1.5*a*(1+s)*(1-s^2) / (1+3*s^2)
で描画する。（ということだ）

s = seq(-3.5, 3.5, by=0.01)
x = 1.5*a*(1-s)*(1-s^2) / (1+3*s^2)
y = 1.5*a*(1+s)*(1-s^2) / (1+3*s^2)
plot(x, y, type="l", asp=1)
abline(-a, -1, col=2, lty=2)
abline(0, 1, col=4, lty=2)
abline(v=0, h=0, col=3, lty=3)

極座標系で描くこともできるが，theta の変閾の設定が若干直感的でない気がした。

theta = seq(-0.6, 2.2, by=0.01)
r = 3 * a * sin(theta) * cos(theta) / (sin(theta)^3 + cos(theta)^3)
x = r*cos(theta)
y = r*sin(theta)
plot(x, y, type = "l", asp=1)
abline(-a, -1, col=2, lty=2)
abline(0, 1, col=4, lty=2)
abline(v=0, h=0, col=3, lty=3)

ベアストウ法による求解

別名ヒチッコック法（というか，ヒチコック・ベアストウ法として覚えていたが）により，実係数の高次代数方程式を解く。虚数解になる場合も含み，n 次式の場合は n 個の解を全て求める。

インターネット上にもプログラム例が見られるが，数値計算的，プログラミング的に見た場合，ちょっとひどいなあというのも（箇所も）見られるので，お手本とはいわないまでも，まずまずのプログラムを示しておこう。

まずは R によるプログラム例，その後に Python によるプログラム例を掲示する。

quadratic.root = function(p, q) {
  d = p ^ 2 - 4 * q
  if (d >= 0) {
    x1 = (-p - ifelse(p >= 0, 1,-1) * sqrt(d)) / 2 # 絶対値の大きいほうの解を求める
    x2 = q / x1 # もう一方は解と係数の関係から得る
  } else {
    x1 = (-p + sqrt(d + 0i)) / 2 # 大抵のプログラミング言語は虚数を取り扱える
    x2 = Conj(x1) # 共役複素数を求める関数もある
  }
  return(c(x1, x2))
}

Bairstow = function(a, EPSILON = 1e-14) {
  a = a / a[1]
  ans = NULL
  repeat {
    n = length(a)
    if (n == 3) { # n = 3, n = 2 の場合の処理は，ループの最初に書かねばならない
      ans = c(ans, quadratic.root(a[2], a[3]))
      break
    } else if (n == 2) {
      ans = c(ans, -a[2] / a[1])
      break
    }
    p = q = 1
    b = c = numeric(n)
    c[1] = b[1] = 1 # a[1], b[1], c[1] は常に 1 である
    repeat {
      b[2] = a[2] - p * b[1]
      c[2] = b[2] - p * c[1]
      for (i in 3:n) {
        b[i] = a[i] - p * b[i - 1] - q * b[i - 2]
        c[i] = b[i] - p * c[i - 1] - q * c[i - 2]
      }
      d = c[n - 2] ^ 2 - c[n - 3] * (c[n - 1] - b[n - 1])
      delta.p = (b[n - 1] * c[n - 2] -  b[n] * c[n - 3]) / d
      delta.q = (b[n] * c[n - 2] - b[n - 1] * (c[n - 1] - b[n - 1])) / d
      if (abs(delta.p) < EPSILON && abs(delta.q) < EPSILON) {
        break
      }
      p = p + delta.p
      q = q + delta.q
    }
    ans = c(ans, quadratic.root(p, q)) # x^2 + px + q = 0 を解く
    a = b[1:(n - 2)] # 元の多項式を x^2 + px + q = 0 で割った多項式の係数
  }
  return(sort(ans))
}

> options(digits = 14)
> Bairstow(c(1.0, -15.0, 85.0, -225.0, 274.0, -120.0))
[1] 1 2 3 4 5
> Bairstow(c(1, 10, 35, 50, 24))
[1] -4 -3 -2 -1
> Bairstow(c(1, 10, 25, 50, 24))
[1] -7.49826796187678+0.0000000000000i
[2] -0.93451222322734-2.0458454872479i
[3] -0.93451222322734+2.0458454872479i
[4] -0.63270759166854+0.0000000000000i
> Bairstow(c(1, 1, 1, 1))
[1] -1+0i  0-1i  0+1i
> Bairstow(c(1, 1, 1))
[1] -0.5-0.86602540378444i -0.5+0.86602540378444i
> Bairstow(c(3, 5))
[1] -1.6666666666667
> Bairstow(c(1,-2, 1))
[1] 1 1
> quadratic.root(-1.000000001, 0.000000001)
[1] 1e+00 1e-09
> quadratic.root(0, 1)
[1] 0+1i 0-1i

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

Python による場合

import scipy as sp

def quadratic_root(p, q):
  d = p**2 - 4 * q
  if d >= 0:
    x1 = (-p - sp.copysign(sp.sqrt(d), p)) / 2
    x2 = q / x1
  else:
    x1 = (-p + sp.sqrt(d + 0j)) / 2
    x2 = x1.conjugate()
  return x1, x2


def Bairstow(a, EPSILON = 1e-14):
  a = sp.ravel(a)
  a = a / a[0]
  ans = []
  while True:
    n = len(a) - 1
    if n == 2:
      ans.extend(quadratic_root(a[1], a[2]))
      break
    elif n == 1:
      ans.append(-a[1] / a[0])
      break
    p = q = 1
    b = sp.ones_like(a)
    c = sp.ones_like(a)
    while True:
      b[1] = a[1] - p * b[0]
      c[1] = b[1] - p * c[0]
      for i in range(2, n + 1):
        b[i] = a[i] - p * b[i - 1] - q * b[i - 2]
        c[i] = b[i] - p * c[i - 1] - q * c[i - 2]
      d = c[n - 2]**2 - c[n - 3] * (c[n - 1] - b[n - 1])
      delta_p = (b[n - 1] * c[n - 2] - b[n] * c[n - 3]) / d
      delta_q = (b[n] * c[n - 2] - b[n - 1] * (c[n - 1] - b[n - 1])) / d
      if abs(delta_p) < EPSILON and abs(delta_q) < EPSILON:
        break
      p += delta_p
      q += delta_q
    ans.extend(quadratic_root(p, q))
    a = b[:n - 1]
  return sp.sort(ans)

>>> Bairstow(sp.array([1, -15, 85, -225, 274, -120.0]))
array([1., 2., 3., 4., 5.])
>>> Bairstow([1, 10, 35, 50, 24])
array([-4., -3., -2., -1.])
>>> Bairstow([1, 10, 25, 50, 24])
array([-7.49826796+0.j        , -0.93451222-2.04584549j,
       -0.93451222+2.04584549j, -0.63270759+0.j        ])
>>> Bairstow([1, 1, 1, 1])
array([-1.+0.j,  0.-1.j,  0.+1.j])
>>> Bairstow([1, 1, 1])
array([-0.5-0.8660254j, -0.5+0.8660254j])
>>> Bairstow([3, 5])
array([-1.66666667])
>>> Bairstow([1, -2, 1])
array([1., 1.])
>>> quadratic_root(-1.000000001, 0.000000001)
(1.0, 1e-09)
>>> quadratic_root(0, 1)
(1j, -1j)

グラフの軸の目盛り設定

中澤さんが

> 出勤前にふと思いついて，trux.Rというコードを書いてみた。これで問題なければ，人口ピラミッド描画パッケージのデフォルト横軸目盛設定アルゴリズムはこれに変えようと思う。

trux <- function(X) {
 SC <- 10^as.integer(log10(X)-1)
 MX <- (X %/% (6*SC))*SC*2
 return(0:4*MX)
}

> trux(132)
[1]   0  40  80 120 160
> trux(495)
[1]   0 160 320 480 640
> trux(8734)
[1]     0  2800  5600  8400 11200
> trux(13044050)
[1] 0.0e+00 4.0e+06 8.0e+06 1.2e+07 1.6e+07
> trux(94822324)
[1] 0.0e+00 3.0e+07 6.0e+07 9.0e+07 1.2e+08

R には pretty という関数が用意されているので，そちらを使うほうがよさそうに思えるのだけど。
コメントは twitter でというようになっているんだけど，私は twitter やんないので。誰か，中澤さんに聞いてみて。

> pretty(c(0, 132))
[1]   0  20  40  60  80 100 120 140
> pretty(c(0, 495))
[1]   0 100 200 300 400 500
> pretty(c(0, 8734))
[1]     0  2000  4000  6000  8000 10000
> pretty(c(0, 13044050))
[1] 0.0e+00 2.0e+06 4.0e+06 6.0e+06 8.0e+06 1.0e+07 1.2e+07 1.4e+07
> pretty(c(0, 94822324))
[1] 0e+00 2e+07 4e+07 6e+07 8e+07 1e+08

ラベルの数は n で指定できるが，必ず指定した数になるわけではない。

> pretty(c(0, 13044050), n=3)
[1] 0.0e+00 5.0e+06 1.0e+07 1.5e+07

その他の引数もある。

ポケモン Go と負の二項分布

よう！ヨッシーじゃないか。歩きながらポケモンやってると危ないよ。

おう！ジュンか。ついつい夢中になっちゃったんだけど，反省しないといけないね。

なんか，おもしろい話はないかい？

そうだねえ。ポケモンといえば，最近ポケモンの写真を何回か撮ると，何回目かに特別なポケモンが写真に写り込んでその特別なポケモンをゲットできるという機能が実装されたんだけど，1 回でゲットできることはほとんどなくて，大抵 2，30 回目にやっとゲットできるようなんだ。ゲットできるまでに何回くらい写真を撮らないといけないかは，1 回あたりのポケモンゲットの確率とどんな関係にあるんだろうね？

「確率 p で起きる事象が k 回起きるまでに，その事象が何回起きなかったか」という分布を負の二項分布というというのがあったね。

負の二項分布か〜。その分布は何で決まるの？ p と k と...

それと，起きなかった回数 x だね。x はトライするたびに変わる。ヨッシーのいっているのは，k = 1 のときだけど，特別なポケモンが写真に写る確率 p は公表されていないよね。

起きなかった回数 x は，経験的に 20 〜 40 回くらいだと思うんだけど，平均的にはどれくらいなんだろうか。

さすがにオレも計算式は暗記していないし，これから講義があるから，それが終わってから昼休みの学食で続きをやろう。統計学の本を持って行くからヨッシーはノートパソコンを持ってこいよ。

じゃあ，後で...

:

さてと。負の二項分布は
f(x) = choose(x+k-1, x) * p^k * q^x
だね。

相変わらず，数式の表現は R の記法を拝借しているんだね（笑）でも，R なら関数が用意されているんじゃないか？

そう。dnbinom という関数がある。dnbinom(x, k, p) で そのときの確率 f(x) が得られる。

x は一つの値でなくて，ベクトルでも良いようだね。x = 0:30, k = 1。p はかなり低い気がするので，p = 0.1 でまずやってみよう。

> (prob = dnbinom(0:30, 1, 0.1))
 [1] 0.100000000 0.090000000 0.081000000 0.072900000 0.065610000
 [6] 0.059049000 0.053144100 0.047829690 0.043046721 0.038742049
[11] 0.034867844 0.031381060 0.028242954 0.025418658 0.022876792
[16] 0.020589113 0.018530202 0.016677182 0.015009464 0.013508517
[21] 0.012157665 0.010941899 0.009847709 0.008862938 0.007976644
[26] 0.007178980 0.006461082 0.005814974 0.005233476 0.004710129
[31] 0.004239116

それじゃよく分からないから，図にして。

じゃあ。

> plot(0:30, prob)


あれ！全然想像と違う。これじゃ，x = 0 のときが一番確率が高い。つまり 1 回目で事象が起きることがもっとも期待されるということだね。経験上は 2，30 回かかると思ってるんだけど。

平均値は k * (1-p) / p，分散は k * (1-p) / p^2 ということだ。平均値は 1.5，分散は 3.75 だ。

p が変わったらどうなるか，やってみよう。

> plot(0:30, dnbinom(0:30, 1, 0.2))


いろいろ変えても同じ傾向だね。x = 0 のときが一番確率が高い。

この分布の形って，k でも変わるよね。k = 5 ぐらいでやってみたら？

> plot(0:30, dnbinom(0:30, 5, 0.2))


あ。こんな感じ。も少し調整してみよう。

これでどうだい。

> plot(0:100, dnbinom(0:100, 8, 0.2))


いい感じ。

平均値は 8 * 0.8 / 0.2 だから 32 だね。これは「事象が発生しなかった回数」だから，何回目に k 回事象が発生したかは 32+8 = 40 回目ということだね。

でも，どうやって解釈すればいいんだろうか？40 回目前後に当たりが出るようにでも設定しているのかなあ？

そうかもしれないけど，ゲームということを考えれば，1 ポイントとるごとにゲットされやすくなっていき，k ポイントとったらゲットできるようになっているのかもしれないね。そう考えれば，この問題は負の二項分布だけで説明できるんじゃないかな？

ポケモンもおもしろいけど，統計学もおもしろいね。

ハハハ。

=== おまけ ===

負の二項分布を R でシミュレーションしてみよう。

sim = function(k, p) {
    count = -k # count は通算何回目に k 回目の事象が起こるかを表している
    repeat {
        count = count + 1
        k = k - (runif(1) <= p) # 事象が発生したら k を 1 減らす
        if (k == 0) { # 0 になったら結果を返す（事象が起きなかった回数）
            return(count)
        }
    }
}

> k = 8
> p = 0.2
> q = 1-p
> z = replicate(10000, sim(k, p))

> mean(z) # シミュレーション結果の平均
[1] 31.9754
> k*q/p # 理論上の平均値
[1] 32

> var(z) # シミュレーション結果の分散
[1] 162.2052
> k*q/p^2 # 理論上の分散
[1] 160

アルキメデスの方法による円周率の近似値

直径 1 の円の外接正 n 角形の周長を Pn，内接正 n 角形の周長を pn とすれば，

n = 3 * 2 ^ k, k = 0, 1, 2, ...

P2n = 2 * Pn * pn

p2n = sqrt(P2n * pn)

ということで，以下のプログラム

import scipy as sp

def Pi_Archimedes():
    print("{0:>2s} {1:>9s} {2:>17s} {3:>17s}".format("k", "n", "Pn", "pn"))
    k = 0 # 正三角形から始める
    Pn = 3 * sp.sqrt(3)
    pn = Pn / 2
    while True:
        n = 3 * 2 ** k
        print("{0:2d} {1:9d} {2:.15f} {3:.15f}".format(k, n, Pn, pn))
        if Pn == pn:
            return Pn
        k += 1
        P2n = 2 * Pn * pn / (Pn + pn)
        p2n = sp.sqrt(P2n * pn)
        Pn, pn = P2n, p2n

Pi_Archimedes()

 k         n                Pn                pn
 0         3 5.196152422706632 2.598076211353316
 1         6 3.464101615137754 3.000000000000000
 2        12 3.215390309173473 3.105828541230249
 3        24 3.159659942097500 3.132628613281238
 4        48 3.146086215131435 3.139350203046867
 5        96 3.142714599645368 3.141031950890509
 6       192 3.141873049979824 3.141452472285462
 7       384 3.141662747056849 3.141557607911857
 8       768 3.141610176604690 3.141583892148318
 9      1536 3.141597034321526 3.141590463228050
10      3072 3.141593748771351 3.141592105999271
11      6144 3.141592927385096 3.141592516692157
12     12288 3.141592722038614 3.141592619365384
13     24576 3.141592670701998 3.141592645033691
14     49152 3.141592657867844 3.141592651450767
15     98304 3.141592654659306 3.141592653055036
16    196608 3.141592653857171 3.141592653456104
17    393216 3.141592653656637 3.141592653556370
18    786432 3.141592653606503 3.141592653581437
19   1572864 3.141592653593970 3.141592653587703
20   3145728 3.141592653590837 3.141592653589270
21   6291456 3.141592653590053 3.141592653589661
22  12582912 3.141592653589857 3.141592653589759
23  25165824 3.141592653589808 3.141592653589784
24  50331648 3.141592653589796 3.141592653589790
25 100663296 3.141592653589793 3.141592653589791
26 201326592 3.141592653589792 3.141592653589792
27 402653184 3.141592653589792 3.141592653589792

3.141592653589792

正 4 億多角形ですと！

 

R によるプログラムもほとんど同じで，以下のように

Pi.Archimedes = function() {
    Pn = 3 * sqrt(3)
    pn = Pn / 2
    while (TRUE) {
        if (Pn == pn) {
            return(Pn)
        }
        tmp = 2 * Pn * pn /  (Pn + pn)
        pn = sqrt(tmp * pn)
        Pn = tmp
    }
}

> options(digits=16)
> Pi.Archimedes()
[1] 3.141592653589792

追記

Julia では

function piarchimedes()
    Pn = 3sqrt(3)
    pn = Pn / 2
    while true
        Pn == pn && return Pn
        tmp = 2Pn * pn /  (Pn + pn)
        pn = sqrt(tmp * pn)
        Pn = tmp
    end
end

piarchimedes() # 3.141592653589792

ロンバーグ積分（その2）

前報の「ロンバーグ積分」には，C で書かれたプログラムを Python に書き換えたものを示した。両言語は共に「0 オリジン」

https://jp.mathworks.com/matlabcentral/fileexchange/34-rombint-m には，MATLAB で書かれたプログラムがある。

# ROMBINT     Numerical evaluation of an integral using the Romberg method.
#
#   Q = rombint(F, A, B) approximates the integral of F(X) from A to B to
#   within a relative error of 1e-10 using Romberg's method of integration.
#   F is a function.  The function must return a vector of output values if
#   a vector of input values is given.
#
#   Q = rombint(F, A, B, DECIMALDIGITS) integrates with accuracy 10^{-DECIMALDIGITS}.
#   ---------------------------------------------------------
#   Author: Martin Kacenak,
#           Department of Applied Informatics and Process Control,
#           Faculty of BERG, Technical University of Kosice,
#           B.Nemcovej 3, 04200 Kosice, Slovak Republic
#   E-mail: ma.kac@post.cz
#   Date:   posted February 2001, updated June 2006.

MATLAB は「1 オリジン」なので，同じく「1 オリジン」である R で書き換えるのが面倒がない。
見かけはちょっと違うが，ほぼ同じ長さのプログラムである。

rombint = function(funfcn, a, b, decdigs = 10) {
  rom = matrix(0, 2, decdigs)
  romall = numeric(2 ^ (decdigs - 1) + 1)
  h = b - a
  romall = funfcn(seq(a, b, by = h / 2 ^ (decdigs - 1)))
  rom[1, 1] = h * (romall[1] + romall[length(romall)]) / 2
  for (i in 2:decdigs) {
    st = 2 ^ (decdigs - i + 1)
    rom[2, 1] = (rom[1, 1] + h * sum(romall[seq(st / 2 + 1, 2 ^ (decdigs - 1), st)])) / 2
    for (k in 1:(i - 1)) {
      rom[2, k + 1] = ((4 ^ k) * rom[2, k] - rom[1, k]) / (4 ^ k - 1)
    }
    rom[1,] = rom[2, ]
    h = h / 2
  }
  return(rom[1, decdigs])
}

Pi = function(x) {
  return(4/(1+x^2))
}

> rombint(Pi, 0, 1)
[1] 3.141592653589793

> pi - rombint(Pi, 0, 1)
[1] 4.440892098500626e-16

ロンバーグ積分

ロンバーグ積分とは，まず合成台形公式で近似し、次に Richardson の補外法で近似をより正確にする。

Romberg's method

に，C で書かれたプログラムがある。

これを，Python に書き換えると以下のようになる。

for をベクトル計算にするなどにより，ずいぶん短くなる？

import scipy as sp

def Romberg(f, a, b, max_steps=1000, acc=1e-14):
    Rp = sp.empty(max_steps)  # previous row
    Rc = sp.empty(max_steps)  # current row
    h = b - a
    Rp[0] = (f(a) + f(b)) * h / 2
    for i in range(1, max_steps):
        h /= 2
        Rc[0] = h * sum(f(sp.arange(1, 2 ** i, 2) * h)) + Rp[0] / 2  # R(i,0)
        for j in range(1, i + 1):
            n_k = 4 ** j
            Rc[j] = (n_k * Rc[j - 1] - Rp[j - 1]) / (n_k - 1)  # R(i,j)
        if i > 1 and abs(Rp[i - 1] - Rc[i]) < acc:
            return Rc[i - 1]
        Rp, Rc = Rc, Rp # swap Rn and Rc
    return Rp[max_steps - 1]  # best guess

def Pi(x):
    return 4 / (1 + x**2)

a = Romberg(Pi, 0, 1, max_steps=9, acc=1e-14)
a　# 3.141592653589794
sp.pi - a  # -8.881784197001252e-16

積分について

積分は，台形公式，シンプソンの公式のほかにブールの公式がある。

目的とする関数は以下のものを例示する

R

Pi = function(x) {
  return(4/(1+x^2))
}

Python

def Pi(x):
    return 4/(1+x**2)

Trapezoidal rule

プログラムは以下のようになる

R

Trapezoidal = function(f, lo, hi, n = 1000) {
  x = seq(lo, hi, length=n+1)
  w = diff(x)[1]
  return(w*(f(x[1])/2 + sum(f(x[2:n])) + f(x[n+1])/2))
}

> pi - Trapezoidal(Pi, 0, 1)
[1] 1.666666662458738e-07

Python

def Trapezoidal(f, lo, hi, n = 1000):
    x = sp.linspace(lo, hi, n+1)
    w = sp.diff(x)[0]
    return w*(f(x[0])/2 + sum(f(x[1:n])) + f(x[n])/2)

sp.pi - Trapezoidal(Pi, 0, 1) # 1.6666666891040904e-07

Simpson's rule

合成シンプソン則のプログラムは以下のようになる

R

Simpson = function(f, lo, hi, n = 1000) {
  x = seq(lo, hi, length=n+1)
  w = diff(x)[1]
  return(w*(f(x[1])+4*sum(f(x[seq(2, n, by = 2)]))+2*sum(f(x[seq(3, n, by = 2)]))+f(x[n+1]))/3)
}

> pi - Simpson(Pi, 0, 1)
[1] 0

Python

def Simpson(f, lo, hi, n = 1000):
    x = sp.linspace(lo, hi, n+1)
    w = sp.diff(x)[0]
    return w*(f(x[0])+4*sum(f(x[1:n:2]))+2*sum(f(x[2:n:2]))+f(x[n]))/3

sp.pi - Simpson(Pi, 0, 1) # -3.1086244689504383e-15

Boole's rule

プログラムは以下のようになる

R

Boole = function (f, lo, hi, n = 1000) {
  n = ceiling(n / 4) * 4
  x = seq(lo, hi, length=n+1)
  w = diff(x)[1]
  return(2*w*(
      7*sum(f(x[seq(1, n, by=4)]))
    +32*sum(f(x[seq(2, n, by=4)]))
    +12*sum(f(x[seq(3, n, by=4)]))
    +32*sum(f(x[seq(4, n, by=4)]))
    + 7*sum(f(x[seq(5, n+1, by=4)])))/45)
}

> pi - Boole(Pi, 0, 1)
[1] 0

Python

def Boole(f, lo, hi, n = 1000):
    n = int(sp.ceil(n / 4) * 4)
    x = sp.linspace(lo, hi, n+1)
    w = sp.diff(x)[0]
    return 2*w*(
        7*sum(f(x[0:n:4]))
        +32*sum(f(x[1:n:4]))
        +12*sum(f(x[2:n:4]))
        +32*sum(f(x[3:n:4]))
        +7*sum(f(x[4:n+1:4])))/45

sp.pi - Boole(Pi, 0, 1) # 0.0

 

ブールの公式はなかなかのものであることがわかった。

 

ルーカス数列

ルーカス数列は，フィボナッチ数列から生成される。

Fibonacci numbers are define by: Fn=Fn-1 + Fn-2

Lucas numbers are define by: Ln=Fn + 2Fn-1

フィボナッチ数列のときと同じようにして，

import scipy as sp

n = 100000
tbl = sp.zeros(10, dtype=int)
a, b = 1, 3
tbl[1], tbl[3] = 1, 1
for i in range(2, n):
  a, b = b, a+b
  tbl[int(str(b)[0])] += 1
print(tbl[1:])

これによって，先頭桁が 1 ~ 9 である項数が以下のようであることが分かる。

[30104 17609 12494  9691  7919  6695  5799  5115  4574]

フィボナッチ数列の場合は

[30103 17610 12494  9690  7918  6695  5798  5117  4575]

だったので，ほとんど同じである。両方ともベンフォードの法則にしたがうといえる。
それにしても，実に良く一致している。