Julia に翻訳--237 モンテカルロ法でeを求める

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

モンテカルロ法でeを求める
https://oku.edu.mie-u.ac.jp/~okumura/python/e-montecarlo.html

ファイル名: montecarlo.jl　　関数名:

翻訳するときに書いたメモ

==========#
#=====
一つの関数にしてJIT化する：
@jit(nopython=True)
def main():
    def ne():
        s = 0
        i = 0
        while s < 1:
            s += random.random()
            i += 1
        return i
    n = 100000000
    t = 0
    for _ in range(n):
        t += ne()
    print(t / n)

%time main()

2.23s になった。
=====#

システムはとても古くて，M1 Mac mini には及びもつかないが

macOS Catalina
vergion 10.15.7
MacBook Pro(Retina, mid 2012)
processor 2.7 GHz quadcore Intel Core i7
mamory 16GB 1600 MHz DDr3

Julia だと,Python とほとんど同じコーディングで

using Statistics

function test()
    function ne()
        s = 0
        i = 0
        while s < 1
            s += rand()
            i += 1
        end
        i
    end
    n = 100000000
    t = 0
    for i = 0:n-1
        t += ne()
    end
    println(mean(t)) # 2.71824778
end

@time test()
#==
2.71822072e8
  3.672514 seconds (36 allocations: 1.109 KiB)
==#

いい勝負ですよ。お金を掛けなくても，Julia に乗り換えるだけで，M1 Mac に太刀打ちできる？？？？

Julia に翻訳--236 十進多倍長計算

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

十進多倍長計算
https://oku.edu.mie-u.ac.jp/~okumura/python/decimal.html

ファイル名: 十進多倍長計算.jl　　関数名:

翻訳するときに書いたメモ

==========#

#> Excelに =0.3-0.2-0.1 と打ち込んでみてください。当然ながら（？）ぴったりゼロになります。

0.3-0.2-0.1　# -2.7755575615628914e-17

#> Excelの計算結果としては，この -2.77556E-17 のほうが正しいのです

#> AppleのNumbersという表計算ソフトの最新版では =74724506/23785549-PI() と打ち込むとぴったり -2.22712210211597E-15 になります

#> Pythonでは decimal というパッケージを使えば任意精度の10進演算ができます。
#> 中略
#> 74724506 / 23785549 - π の値を求めてみましょう：
#> 中略
#> 結果が -2.22712210211597472...E-15 のように出力されます。

Julia では

74724506/23785549-π　# -2.220446049250313e-15　？？

以下のようにすれば
big(74724506) / big(23785549) - π # -2.227122102115974720471985452108342639006776121601087390928614849464470859988176e-15

#> おまけ）1 / 0.9899 を計算して，小数点以下を2桁ずつ区切って読んでみてください。フィボナッチ数列になっています
（10000/9899 を計算するとなぜ商 (1.0102030508) にフィボナッチ数列が現れるのですか？）\https://jp.quora.com/10000-9899-%E3%82%92%E8%A8%88%E7%AE%97%E3%81%99%E3%82%8B%E3%81%A8%E3%81%AA%E3%81%9C%E5%95%86-1-0102030508-%E3%81%AB%E3%83%95%E3%82%A3%E3%83%9C%E3%83%8A%E3%83%83%E3%83%81%E6%95%B0%E5%88%97%E3%81%8C%E7%8F%BE%E3%82%8C

1 / big"0.9899" # 1.01020305081321345590463683200323264976260228305889483786241034447924032730578

1.01 02 03 05 08 13 21 34 55 90463683200323264976260228305889483786241034447924032730578

Julia に翻訳--235 PISA「盗難事件」問題のシミュレーション 二項分布

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

PISA「盗難事件」問題のシミュレーション
https://oku.edu.mie-u.ac.jp/~okumura/python/pisa-robbr.html

ファイル名: PISA.jl　　関数名:

翻訳するときに書いたメモ

==========#

#> 全部で1024件の盗難事件が1998年・1999年の2年間に起こった場合，どちらの年も盗難の起きやすさは同じと仮定して（帰無仮説），両者の件数が8件以上異なる確率（p 値）をシミュレーションで求める。

#> import numpy as np

#> rng = np.random.default_rng()
#> a = rng.integers(0, 2, size=1024)

#> 両年の件数の差は 2 * abs(sum(a) - 512) である。これが8以上なら True，8未満なら False を返す関数を作ってみる。

#> def f():
#>     a = rng.integers(0, 2, size=1024)
#>     return 2 * abs(sum(a) - 512) >= 8

#> これの返す値を要素とする長さ10000の配列を作り，True の割合を調べる。

#> b = [f() for _ in range(10000)]
#> sum(b) / 10000

なにか，難しそう。

#> rng = np.random.default_rng()
#> a = rng.integers(0, 2, size=1024)

って？

要するに，ランダムな 0/1 の二値データを 1024 個発生させる？0/1 の差は 2 * abs(sum(a) - 512)

1024 個のランダムな 0/1 二値データ（0 も 1 も同じ確率である）は，

a = rand(0:1, 1024)

で生成される。

1 であるのは

sum(a) # sum(a .== 1) でもあるが

0 であるのは

1024 - sum(a) # sum(a .== 0)

この差は

sum(a) - (1024 - sum(a))

すなわち

2 * sum(a) - 1024

絶対値を取ると

abs(2 * sum(a) - 1024)

すなわち

2 * abs(sum(a) - 512)

> これが8以上なら True，8未満なら False を返す関数

というのも，なぜ 8 以上か未満かというのもいきなりなので，まずは 2 * abs(sum(a) - 512) の分布を調べてみよう。

1024 個のランダムな 0/1 データで，0 と 1 の個数の差の絶対値は，

f() = 2 * abs(sum(rand(0:1, 1024)) - 512) # a = rand(0:1, 1024)

これを 10000 回調べてみよう。

using Random; Random.seed!(123)
results = [f() for i = 1:10000];

results の分布は，たとえば

maxx = maximum(results)

0〜maxxの範囲にある（偶数しかとらないが）。

集計すると

frequency = zeros(maxx+1);
[frequency[i+1] += 1 for i in results];

using Plots
bar(frequency, label="")
vline!([8], label="")
savefig("fig.png")

まあ，差がぴったり0というのはさすがにそう多くはないが差が 20 以内というのはけっこうある。
では，差が 8 以上というのがどれくらいあるかというと

]frequency[1:8]
sum(frequency[1:8]) # 1693.0
sum(results .<= 6)  # 1693
sum(results .>= 8)  # 8307

今の場合は 8307 回であった。つまり 10000 回のうちの 83.07% ということで，「差が 8 以上」ということは，よくあることということになる。

#> この値を p 値という。