Julia でアニメーション GIF，モンテカルロ法で円周率を求める

R には動画画像ファイルの作成のために animation ライブラリがあるが，Julia にも動画 gif ファイルを作成する機能がある。

試しに，モンテカルロ法で π を求める動画を作ってみた。

以下は，2つの gif ファイルを左右に並べて表示しているので（環境によっては上下になっているかも知れない）双方で同期が取れていないのだが，これを一つのファイルにする方法がわからない。色々やってみたのだけど，どれも旨くいかなかった。やりかたをご存じの方お教えくださればうれしいです。

using Random, Plots

function MonteCarlo(n = 100)
    Random.seed!(777)
    x = range(0, 1, length=1000)
    y = @. sqrt(1-x^2);
    xy = rand(n, 2);
    r2 = vec(sum(xy .^ 2, dims=2));
    inout = (r2 .< 1) .+ 0;
    cumpi = cumsum(inout) * 4 ./ collect(1:n);
    pyplot(grid=false, label="", size=(300, 300))
    plt = plot(x, y, color=:black, xlims=(0, 1.01), ylims=(0, 1), aspect_ratio=1)
    plt = plot!([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], color=:black)
    anim = @animate for i = 1:div(n, 4) * 4
        plt = scatter!([xy[i, 1]], [xy[i, 2]],
            color= inout[i] == 1 ? :red : :blue, markerstrokewidth=0)
    end
    gif(anim, "MonteCarlo.gif", fps=10)
    yrange = extrema(cumpi).+ (-0.1, 0.1)
    plt2 = plot([0, n], [π, π], color=:black, ylims=yrange,
        xlabel="number of trial", ylabel="approximated value")
    #plt2 = plot!([0, 0], [yrange[1], yrange[2]], color=:black)
    #plt2 = vline!([0], color=:black)
    plt2 = hline!([π], color=:blue)
    anim2 = @animate for i = 2:div(n, 4) * 4
        plt2 = plot!([i-1, i], [cumpi[i-1], cumpi[i]], color=:black)
    end
    gif(anim2, "MonteCarlo2.gif", fps=10)
end

MonteCarlo(300)

 

自由落下の危険性

自由落下における速度，距離，時間

重力加速度(g) ≒ 9.8 m/s^2

時間(t;s)と速度(v;m/s)

次元解析　(m/s) = (m/s^2) * （s)

v = gt

時間(t;s)と落下距離(h;m)

次元解析　(m) = (m/s^2) * （s^2)

h = gt^2 / 2

落下距離(h;m)に必要な時間(t;s)

次元解析　（s) = √((m) / (m/s^2))

t = √(2h / g)

落下距離(h;m)のときの速度(v;m/s)

次元解析　(m/s) = √(m/s^2 * (m))

v = √(2gh)

単位を km/h にするには

v = √(2gh) * 60^2 / 1000

using Plots
pyplot(size=(300, 200), grid=true, label="")
g = 9.8
vt(t) = g*t
h(t) = g * t^2/2
t(h) =sqrt(2h/g)
v(h) = sqrt(2g*h)*3.6
height = 0:0.1:30
velocity = round(v(14), digits=1)
plot(height, v.(height), tick_direction=:out)
annotate!(12, 75, text("高さ = 14m, 速度 = $(velocity)km/h", 8))

高さ14mから落ちるのは，時速60キロで衝突する衝撃とほぼ同じということだ。

R で ペル方程式，Pell's equation

Mac OS Catalina,  バージョン 10.15.7 で gmp 0.6−2 が動かない。爆弾（エラー）が出る。

いろいろ探して，

M1チップ搭載のMacだとRのgmpパッケージが動かない？
https://rion778.hatenablog.com/entry/2021/03/02/214104

に，

install.packages("gmp", type = "source")

とすればよいと書いてあった。
私のマシンはポンコツなので，M1 なんかでないけど，とにかくやってみたら，動くようになった。

やってみたかったのは，R での Pell's equation プログラムを確かめることだ。

Pell's equation
https://rosettacode.org/wiki/Pell%27s_equation

ここには，いろいろな言語で書かれた，ペル方程式の解を求めるプログラムがのっているが，R によるプログラムがないので，以下のように翻訳してみた。
a, b = b, c みたいな構文がないのでちょっとダサいプログラムになる。

library(gmp)
pell = function(n) {
    x = as.bigz(floor(sqrt(n)))
    y = x
    z = 1
    r = 2 * x
    e1 = 1
    e2 = 0
    f1 = 0
    f2 = 1
    while(TRUE) {
        y = r * z - y
        z = (n - y * y) %/% z
        r = (x + y) %/% z
        temp = e1
        e1 = e2
        e2 = e2 * r + temp
        temp = f1
        f1 = f2
        f2 = f2 * r + temp
        a = f2
        b = e2
        temp = b
        b = a
        a = a * x + temp
        if (a * a - n * b * b == 1) {
            return(c(a, b))
        }
    }
}

pell(61)
# Big Integer ('bigz') object of length 2:
# [1] 1766319049 226153980

as.bigz("1766319049")^2 - 61 * as.bigz("226153980")^2
# 1

pell(109)
# Big Integer ('bigz') object of length 2:
# [1] 158070671986249 15140424455100

pell(181)
# Big Integer ('bigz') object of length 2:
# [1] 2469645423824185801 183567298683461940

pell(277)
# Big Integer ('bigz') object of length 2:
# [1] 159150073798980475849 9562401173878027020

無事に動いているようだ。

図形問題の解（その 2）

図形問題の解
https://blog.goo.ne.jp/r-de-r/e/6bd08cdd43e3a42155ce0ac444b414f0

鶏を割くに牛刀を以てす　ではあるが，

件の図で，
f を原点(0, 0) として，
fb を f = 5/2 *x
fd を g = x / 4
bd を h = -2x + 9
とし，f - g, h - g の定積分を求め，2倍する

Julia の SymPy でやってみよう

using SymPy

@syms f g h x

f = 5/2 * x
g = x / 4
h = -2x + 9
2(integrate(f - g, (x, 0, 2)) + integrate(h - g, (x, 2, 4)))　# 18

図形問題の解

https://www.youtube.com/watch?v=51Ru7F2uMms&feature=emb_rel_end

小学生でも解けるかな？　っていうか，私が小学生並みなだけ。

下図のように頂点の記号を定める。平行四辺形の2つの頂点 b, d を結ぶ。ab の延長線とed の延長線の交点が c。af, ce は垂直線，ac, ef は水平線。

三角形abfの面積S1は5。三角形bcdの面積S2は4。三角形defの面積S3は2。四角形acef の面積は20。

三角形bdfの面積 = 20 - 5  - 4 - 2 = 9

S = 三角形bdfの面積  × 2 = 9 × 2 = 18

R の bincombinations() を Julia に移植

R の bincombinations()

R の bincombinations(p) は，0/1 からなる 2^p 個の長さ p のベクトルを生成する。

Julia に移植すると以下のようになる。

function bincombinations(p)
    ary = Array{Int,2}(undef, 2^p, p);
    for n = 1:p
        ary[:, n] = repeat(vcat(fill(0, 2^(p - n)), fill(1, 2^(p - n))), 2^(n-1))
    end
    ary
end

p = 3 のときは，以下のような 2^p 行 p 列の行列を返す。

bincombinations(3)

8×3 Matrix&#123;Int64&#125;:
 0  0  0
 0  0  1
 0  1  0
 0  1  1
 1  0  0
 1  0  1
 1  1  0
 1  1  1

行ベクトルは，長さが p で，全ての 0/1 の出現パターンを返す。

例えば，長さが p の場合に，各行の 0 の個数が偶数であるのは何通りあるかというような問題を解くときに利用することができる。

function f(p)
    ary = bincombinations(p)
    sum(sum(ary, dims=2) .% 2 .== 0)
end

f(5) # 16

Julia では，`Iterators.product` を使うのが普通かも知れない。

1 つの方法は，`Iterators.product` を使うが，実際にメモリ上に上述のような行列を作る。

function f2(p)
    all_perm(x, n) = vec(map(collect, Iterators.product(ntuple(_ -> x, n)...)))
    count = 0
    for i in all_perm([0, 1], p)
        count += sum(i) % 2 == 0
    end
    count
end

f2(5) # 16

もう 1 つの方法では，実際にメモリ上に行列を作らない。

function f3(p)
    all_perm(x, n) = Iterators.product([x for i = 1:n]...)
    count = 0
    for i in all_perm([0, 1], p)
        count += sum(i) % 2 == 0
    end
    count
end

f3(5) # 16

いずれも p = 20 程度までは，数秒で答えが得られる。
しかし，p が 23 にもなると，`f2()` と `f3()` は相当な時間を要することになる。

もともとは，メモリー上に作らないのでメモリー節約ということであっただろうが，1要素（1行）ごとに処理しているのもあって時間が掛かるのであろう。
どうせ，p はそんなに大きくないのだから，メモリー上に行列を前もって作る `f()` のほうが優れていると結論できるであろう。

p を 5〜25 として，各関数ごとに 10 回ずつ処理時間を計測し，その中央値を図に描くと以下のようになった。

CPUtime = [
    5 5.00E-06 1.30E-05 1.60E-05
    6 7.00E-06 2.00E-05 2.85E-05
    7 9.00E-06 3.60E-05 4.50E-05
    8 1.65E-05 6.85E-05 0.0001225
    9 4.35E-05 0.0001375 0.0003095
    10 8.55E-05 0.0002665 0.0006935
    11 0.000143 0.0005445 0.0016065
    12 0.0002985 0.0011415 0.0037025
    13 0.0009345 0.00226  0.0032525
    14 0.0020035 0.00456  0.005561
    15 0.004387 0.00938  0.0148675
    16 0.009444 0.0702565 0.087456
    17 0.0192205 0.1532295 0.18541
    18 0.0369055 0.364271 0.372975
    19 0.0873515 0.757955 0.7971625
    20 0.1958075 1.6105315 1.5161105
    21 0.654857 3.273466 3.0952965
    22 2.5566305 8.1765115 7.4054555
    23 6.919126 18.127678 22.38751
    24 14.1758585 44.117498 81.814251
    25 42.42869 136.028486 147.6340555]

using Plots
pyplot(size=(300, 300), guidefontsize=6, label="")
plt = plot(CPUtime[:, 1], CPUtime[:, 2:4], grid=false,
    tick_direction=:out, xlims=(4.5, 27), # yscale=:log10,
    markershape=:circle, markersize=4, markerstrokewidth=0,
    xlabel="p", ylabel="CPU time(seconds)")
annotate!.([25.5], CPUtime[21, 2:4], text.(["f()", "f2()", "f3()"], 8, :left));

plt2 = plot(CPUtime[:, 1], CPUtime[:, 2:4], grid=false,
    tick_direction=:out, xlims=(4.5, 27), yscale=:log10,
    markershape=:circle, markersize=4, markerstrokewidth=0,
    xlabel="p", ylabel="CPU time(seconds)")
annotate!.([25.5], CPUtime[21, 2:4], text.(["f()", "f2()", "f3()"], 8, :left));
plot(plt, plt2, size=(500, 250))

自明な？クイズを SymPy で解くと？？な結果に

https://blog.goo.ne.jp/esp_spicy_sp/1?st=0

で，要するに (x - a)(x - b)…(x - z) を展開するとどうなるかということで，答えは 0，なぜならば途中に (x-x) すなわち 0 があるでしょうということで，まあ，トンチ問題，引っかけ問題ではあるのだけど。

しかし，これを Julia の SymPy で解こうとすると，「え！！」となるかもしれないのだった？

Julia の場合まあ，途中省略で。

using SymPy

@syms a b c d e f g h i j k l m n o p q r s t u v w x y z

expand((x-a)(x-x))  # -a 何だと！！
expand((x-a)*(x-x)) # 0 だよね

Python の場合

from sympy import *
var('x a b c')
expand((x-a)(x-x))  # エラーになる。すなわち TypeError: 'Add' object is not callable
expand((x-a)*(x-x)) # 0 だよね

Julia では，

expand((x-1)(x-5))  # が x - 6 になるんだよ...

expand((x-1)*(x-5)) # は，ちゃんと，x^2 - 6x + 5 だよ

Julia で ペル方程式，Pell's equation

前報の無理数の連分数近似というのは，それ自体はあまりありがたいものでもない。

しかし，ペル方程式 Pell's equation を解くときに力を発揮する。

ペル方程式とは，

n を平方数ではない自然数として、未知整数 x, y について

    x^2 − n*y^2 = 1

の形のディオファントス方程式(整係数多変数高次不定方程式)である。

1 つの解 (x, y) が得られたとき，

    x_k + y_k √n = (x + y √n)^k

は，全てペル方程式の解になる。

最小解は、√n の連分数展開の結果を利用する方法がある。
周期が m のとき，a_0, a_1, ..., a_m とすると，a_m までの近似分数を P/Q としたとき，P, Q がペル方程式の解である。
ただし，m が奇数の場合は得られるのは x^2 − n*y^2 = −1 の解なので，それを 2 乗する。または，二倍の周期，すなわち a_2m までの近似分数 P/Q の P, Q が解である。

function regularcontinuedfractions2(d, n=10) # n は a0 も含めた項数
    x = sqrt(d)
    v = 1
    u = a = floor(Int, x)
    an = [a]
    println("v = $v, a = $a, u = $u")
    for i = 1:n-1
        v2 = (d - u^2) / v
        a2 = floor(Int, (x + u) / v2)
        u2 = a2 * v2 - u
        v, a, u = v2, a2, u2
        println("v = $v, a = $a, u = $u")
        append!(an, a)
    end
    an
end

function approximation2(a)
    p0 = BigInt(a[1])
    q0 = BigInt(1)
    println("0: $p0 / $q0 = $(p0 / q0)")
    p1 = BigInt(a[1]*a[2] + 1)
    q1 = BigInt(a[2])
    println("1: $p1 / $q1 = $(p1 / q1)")
    result = [(p0, q0), (p1, q1)]
    for i = 3:length(a)
        p2 = BigInt(a[i]) * p1 + p0
        q2 = BigInt(a[i]) * q1 + q0
        println("$(i-1): $p2 / $q2 = $(p2/q2)")
        append!(result, [(p2, q2)])
        p0, q0, p1, q1 = p1, q1, p2, q2
    end
    result
end

a = regularcontinuedfractions2(61, 25);
for i = 1:length(a)
    println("a[$(i-1)] = $(a[(i)])")
end
#=
a[0] = 7
a[1] = 1
a[2] = 4
a[3] = 3
a[4] = 1
a[5] = 2
a[6] = 2
a[7] = 1
a[8] = 3
a[9] = 4
a[10] = 1
a[11] = 14
a[12] = 1
a[13] = 4
a[14] = 3
a[15] = 1
a[16] = 2
a[17] = 2
a[18] = 1
a[19] = 3
a[20] = 4
a[21] = 1
a[22] = 14
a[23] = 1
a[24] = 4
=#

周期は 11（1,4,3,1,2,2,1,3,4,1,14） で，奇数なので，a[22] までの解を用いる。

b = approximation2(a);
for i = 1:25
    println("$i, $(b[i][1]), $(b[i][2]), $((b[i][1])^2 - (b[i][2])^2 * 61)")
end
#=
1, 7, 1, -12
2, 8, 1, 3
3, 39, 5, -4
4, 125, 16, 9
5, 164, 21, -5
6, 453, 58, 5
7, 1070, 137, -9
8, 1523, 195, 4
9, 5639, 722, -3
10, 24079, 3083, 12
11, 29718, 3805, -1
12, 440131, 56353, 12
13, 469849, 60158, -3
14, 2319527, 296985, 4
15, 7428430, 951113, -9
16, 9747957, 1248098, 5
17, 26924344, 3447309, -5
18, 63596645, 8142716, 9
19, 90520989, 11590025, -4
20, 335159612, 42912791, 3
21, 1431159437, 183241189, -12
22, 1766319049, 226153980, 1
23, 26159626123, 3349396909, -12
24, 27925945172, 3575550889, 3
25, 137863406811, 17651600465, -4
=#

a[22] までの連分数近似では，(P, Q) = (x, y) = (1766319049, 226153980)

1766319049^2 - 61 * 226153980^2 # 1

「x_k + y_k √n = (x + y √n)^k」というのは，単純な数値計算をするのではない。

(x + √n * y)^2 # 1.2479531931441058e19 という1つの数になる 

(x + √n * y)^2 を数式として展開すると
x^2 + n*y^2 + 2xy√n
で，新たな解 x は x^2 + n*y^2，y は 2xy である。

つまり，√n を含まない項（の和）が x，√n を含む項の係数が y ということである。

x = 1766319049^2 + 61 * 226153980^2 # 6239765965720528801
y = 2 * 1766319049 * 226153980      # 798920165762330040

つまり，(x, y) = (6239765965720528801, 798920165762330040) である。

x^2 - 61 * y^2 # 1

Julia で連分数展開，連分数近似

連分数展開とは，たとえば x の平方根を以下のような連続する分数の形式で表現することである。

√x = a0 + 1/(a1 + 1/(a2 + 1/(a3 + 1/(a4 + 1 / (a5 + ε)))))
εはさらに 1/(a6 + 1/(a7 + 1/...)) と同じように無限に繰り返されるが，だんだんと近似に寄与する割合が小さくなるので，いつかの時点で ε = 0 とすることで近似できる。

例えば x = 31 で a0 から a8 までを使うと，以下のような連分数による近似が得られる。
√31 ≒ 5 + 1/(1 + 1/(1 + 1/(3 + 1/(5 + 1 / (3 + 1/(1 + 1/(1 + 1/(10 + 0))))))))
    = 16063 / 2885 = 5.5677642980935875

なお，全ての分数の分子が 1 であるような場合は，正則連分数とよぶ。

√31 の連分数展開は以下のように行われる。

d = 31
z = sqrt(d) # 5.5677643628300215

1. z の整数部を取る。これが a0 である
    a0 = 5
2. z の小数部の逆数（0.5677643628300215）をとり，新たな z とする
    z = 1 / 0.5677643628300215 # 1.7612940604716716
3. z の整数部を取る。これが a1 である
    a1 = 1
4. z の小数部（0.7612940604716716）の逆数をとり，新たな z とする
    z = 1 / 0.7612940604716716 # 1.3135528725660022
この後も 3.，4. の手順を繰り返す（結局は 1. 2. でもあるが）をくりかえす。

なお，計算の途中で z が 2. のときの z の小数部と等しくなることがある。
そうすると，それ以降の z の整数部は以前のz の整数部の繰り返しになることになる。
√31 の展開の場合は，5, 1, 1, 3, 5, 3, 1, 1, 10, 1, 1, 3, 5, 3, 1, 1, 10, 1, 1, ...
のようになり，1, 1, 3, 5, 3, 1, 1, 10 が循環する。

手順の中に 「z の小数部」が必要になるが，手順を繰り返していくと精度に問題が生じる可能性がある。
そこで，対処法の一つは長精度実数（Julia では BigFloat) を使うことである（問題がないこともあるし，不十分なこともある）。

この計算プログラムでは，平方根の連分数近似だけではなく，引数に与えられる数値（有理数，無理数とも）の連分数近似を行うことができる。

function regularcontinuedfractions(z, n=10) # n は a0 も含めた項数
    an = []
    for i =0:n-1
        ai = Int(floor(z))
        append!(an, ai)
        println("$i: a[$i] = $ai,  z - ai = $(z - ai)")
        z = 1 / (z - ai)
    end
    an
end

a = regularcontinuedfractions(sqrt(big(31)), 11)
#=
0: a[0] = 5,  z - ai = 0.5677643628300219221194712989185495204763933775704143039684325856035898392542376
1: a[1] = 1,  z - ai = 0.7612940604716703203532452164864249200793988962617357173280720976005983065423672
2: a[2] = 1,  z - ai = 0.3135528725660043844238942597837099040952786755140828607936865171207179678508544
3: a[3] = 3,  z - ai = 0.1892547876100073073731570996395165068254644591901381013228108618678632797513434
4: a[4] = 5,  z - ai = 0.2838821814150109610597356494592747602381966887852071519842162928017949196290188
5: a[5] = 3,  z - ai = 0.5225881209433406407064904329728498401587977925234714346561441952011966130611749
6: a[6] = 1,  z - ai = 0.9135528725660043844238942597837099040952786755140828607936865171207179679371505
7: a[7] = 1,  z - ai = 0.09462739380500365368657854981975825341273222959506905066140543093393163977229678
8: a[8] = 10,  z - ai = 0.5677643628300219221194712989185495204763933775704143039684325856035898508027295
9: a[9] = 1,  z - ai = 0.7612940604716703203532452164864249200793988962617357173280720976005982707171392
10: a[10] = 1,  z - ai = 0.3135528725660043844238942597837099040952786755140828607936865171207180296644554

11-element Vector{Any}:
  5
  1
  1
  3
  5
  3
  1
  1
 10
  1
  1
=#

もう一つの対処法は，z の小数部を使わない方法である。
for 文の中の 2 行目，floor(Int, (x + u) / v2) で整数部をとるが，小数部はこのプログラムのどこにも出てこない。
ただし，以下のプログラムは，引数 d の平方根 √d の連分数近似に限定される。

function regularcontinuedfractions2(d, n=10) # n は a0 も含めた項数
    x = sqrt(d)
    v = 1
    u = a = floor(Int, x)
    an = [a]
    println("v = $v, a = $a, u = $u")
    for i = 1:n-1
        v2 = (d - u^2) / v
        a2 = floor(Int, (x + u) / v2)
        u2 = a2 * v2 - u
        v, a, u = v2, a2, u2
        println("v = $v, a = $a, u = $u")
        append!(an, a)
    end
    an
end

rcf = regularcontinuedfractions2(31)
#=
v = 1, a = 5, u = 5
v = 6.0, a = 1, u = 1.0
v = 5.0, a = 1, u = 4.0
v = 3.0, a = 3, u = 5.0
v = 2.0, a = 5, u = 5.0
v = 3.0, a = 3, u = 4.0
v = 5.0, a = 1, u = 1.0
v = 6.0, a = 1, u = 5.0
v = 1.0, a = 10, u = 5.0
v = 6.0, a = 1, u = 1.0

10-element Vector{Int64}:
  5
  1
  1
  3
  5
  3
  1
  1
 10
  1
=#

連分数近似の結果

次に，連分数展開の結果による近似値を計算する関数を定義しよう。
regularcontinuedfractions() が返す結果ベクトル a_i, i = 0, 1, ..., n-1 を与える。

function approximation(a)
    x = a[end]
    for i in reverse(a[2:end-1])
        x = i + 1/x
    end
    a[1] + 1 / x
end

rcf = regularcontinuedfractions2(31, 9) # n は a0 も含めた項数
#=
v = 1, a = 5, u = 5
v = 6.0, a = 1, u = 1.0
v = 5.0, a = 1, u = 4.0
v = 3.0, a = 3, u = 5.0
v = 2.0, a = 5, u = 5.0
v = 3.0, a = 3, u = 4.0
v = 5.0, a = 1, u = 1.0
v = 6.0, a = 1, u = 5.0
v = 1.0, a = 10, u = 5.0

9-element Vector{Int64}:
  5
  1
  1
  3
  5
  3
  1
  1
 10
=#

approximation(rcf) # 5.5677642980935875

途中の近似値を計算するための分数を全て返す関数を定義しよう。
連分数近似 regularcontinuedfractions() が返す結果ベクトル a_i, i = 0, 1, ..., n-1 を与える。
(p(i), q(i)), i = 0, 1, ... のタプルのベクトル。
i までの結果を使う近似値は p(i) / q(i) で得られる。

function approximation2(a)
    p0 = BigInt(a[1])
    q0 = BigInt(1)
    println("0: $p0 / $q0 = $(p0 / q0)")
    p1 = BigInt(a[1]*a[2] + 1)
    q1 = BigInt(a[2])
    println("1: $p1 / $q1 = $(p1 / q1)")
    result = [(p0, q0), (p1, q1)]
    for i = 3:length(a)
        p2 = BigInt(a[i]) * p1 + p0
        q2 = BigInt(a[i]) * q1 + q0
        println("$(i-1): $p2 / $q2 = $(p2/q2)")
        append!(result, [(p2, q2)])
        p0, q0, p1, q1 = p1, q1, p2, q2
    end
    result
end

rcf = regularcontinuedfractions2(31, 9) # n は a0 も含めた項数
b = approximation2(rcf)
#=
0: 5 / 1 = 5.0
1: 6 / 1 = 6.0
2: 11 / 2 = 5.5
3: 39 / 7 = 5.571428571428571428571428571428571428571428571428571428571428571428571428571419
4: 206 / 37 = 5.567567567567567567567567567567567567567567567567567567567567567567567567567558
5: 657 / 118 = 5.567796610169491525423728813559322033898305084745762711864406779661016949152556
6: 863 / 155 = 5.567741935483870967741935483870967741935483870967741935483870967741935483870972
7: 1520 / 273 = 5.567765567765567765567765567765567765567765567765567765567765567765567765567756
8: 16063 / 2885 = 5.567764298093587521663778162911611785095320623916811091854419410745233968804128

9-element Vector{Tuple{BigInt, BigInt}}:
 (5, 1)
 (6, 1)
 (11, 2)
 (39, 7)
 (206, 37)
 (657, 118)
 (863, 155)
 (1520, 273)
 (16063, 2885)
 =#

b[7][1] / b[7][2] # 5.567741935483870967741935483870967741935483870967741935483870967741935483870972

regularcontinuedfractions() は無理数も連分数近似できる。

π の近似

a = regularcontinuedfractions(π)
approximation2(a)
big(π)
#=                    3.141592653589793238462643383279502884197169399375105820974944592307816406286198
9: 1146408 / 364913 = 
 3.14159265359140397848254241421927966391989323482583519907484797746312134673195
=#

ネイピア数

a = regularcontinuedfractions(ℯ, 50);
approximation2(a)
big(ℯ)
#=
 2.718281828459045235360287471352662497757247093699959574966967627724076630353555
9: 1457 / 536 =
 2.718283582089552238805970149253731343283582089552238805970149253731343283582103
49: 3471379623627236787997 / 1277049196033919722542 =
 2.718281828459045067661362633906433282499985585104088805075501397944665141326474
=#

有理数の連分数近似

a = regularcontinuedfractions(big"1.2345", 10);
approximation(a) # 1.2345
approximation2(a)
# 9: 2469 / 2000 = 1.234499999999999999999999999999999999999999999999999999999999999999999999999991

枝葉末節

qiita.com/tags/python で，i += 1 より i = i+1 のほうが速いというのがあったが，そもそも，10000000 回！！もやって所要時間が 0.69685sec と 0.67206sec なんだから，誤差範囲でしょう。

ちなみに，Julia でやったら，誤差範囲で同じようなものだった。
コンパイルするから，記述の差は吸収されてしまうのだろうなあ。

@time begin
        i = 0
        while i < 10_000_000
          i += 1
        end
      end
#   0.520408, 0.507725, 0.509720, 0.508239, 0.530775 などなど

@time begin
        i = 0
        while i < 10_000_000
          i = i + 1
        end
      end
#   0.503195, 0.521761, 0.483759, 0.513919, 0.516283　などなど

繰り返すが，こんな差は，鼻くそみたいなものだ。（^_^;）

クヌースは（C 言語ではあるが） += を推奨している。

abcdefghijklmn += 1

と

abcdefghijklmn2 = abcdefghijklmn + 1

を見極める余計な注意力を必要としないことをよしとするのだ。

 

Python での実行例

Python では i = i + 1 のほうが速いだろうという考察もあるようだが。

from __future__ import print_function
import timeit

print("i+=1\t", timeit.timeit("while i<10000000: i+=1", setup="i=0"))
print(" i=i+1\t", timeit.timeit("while i<10000000: i=i+1", setup="i=0"))

print("a[0]+=1\t", timeit.timeit("while a[0]<10000000: a[0]+=1", setup="a=[0]"))
print("a[0]=a[0]+1\t", timeit.timeit("while a[0]<10000000: a[0]=a[0]+1", setup="a=[0]"))

#=
i+=1  0.6968548460000079
i=i+1  0.6720683110000039
a[0]+=1  1.513158634000007
a[0]=a[0]+1  1.4887571600000058
=#

Julia の小ネタ--027 素数に関連する関数たち

Julia の Primes モジュールは，R や Python で gmp パッケージでサポートされているものと同じだ。

using Primes

素因数分解　factor(n::Integer), factor(ContainerType, n::Integer)

結果の表示法

特に指定しない場合は，数式として見やすい結果が得られる。
540 = 2^2 * 3^3 * 5

factor(540)
# 2^2 * 3^3 * 5

ベクトルとして結果を得るためには ContainerType を Vector とする。同じ素因数が複数個ある場合は，複数個数分の要素として含まれる。

factor(Vector, 540)
#=
6-element Vector{Int64}:
 2
 2
 3
 3
 3
 5
=#

集合として結果を得る場合は，ContainerType を Set とする。この場合は，ユニークな素因数が列挙されるだけであり，それぞれの素因数が何個あったかの情報は失われる。

factor(Set, 540)
#=
Set{Int64} with 3 elements:
  5
  2
  3
=#

キーをソートした辞書型として結果を得るる場合は，ContainerType を DataStructures.SortedDict とする。

factor(DataStructures.SortedDict, 540)
#=
SortedDict{Int64, Int64, Base.Order.ForwardOrdering} with 3 entries:
  2 => 2
  3 => 3
  5 => 1
=#

大きな数の素因数分解

n = 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 29 * 31 * 37 # 7420738134810
factor(n) # 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 29 * 31 * 37

factor(100000000000000000000000000007) # 53 * 827 * 2281490269444000821336497

factor(2281490269444000821336497) # 2281490269444000821336497
factor(big"2281490269444000821336497") # 2281490269444000821336497

factor() の逆関数　prodfactors()

factor() の結果を全部掛け合わすと，元の数になる。

n = 540
prod(factor(n)) # 540
prodfactors(factor(Vector, n)) # 540
prodfactors(factor(DataStructures.SortedDict, n)) # 540
prodfactors(factor(Set, n)) # 30 これは，正しい答えではない

範囲内の素数生成　primes([lo,] hi)

primes(20, 50)
#=
7-element Vector{Int64}:
 23
 29
 31
 37
 41
 43
 47
=#

次の素数　nextprime(n::Integer, i::Integer=1; interval::Integer=1)

n より小さくない i 番目に小さい素数を返す。言い回しがまだるっこしいが，指定された数が素数だった場合，次の（1番目に小さい）素数はその数自体であるということ。

nextprime(100000000000000000000000000000) # 100000000000000000000000000319
factor(100000000000000000000000000319)    # 100000000000000000000000000319

nextprime(100000000000000000000000000319) # 100000000000000000000000000319

前の素数　prevprime(n::Integer, i::Integer=1; interval::Integer=1)

n より大きくない i 番目に大きい素数

2 つ先の素数

nextprime(100000000000000000000000000000, 2) # 100000000000000000000000000379

1 つ前の素数

prevprime(100000000000000000000000000378)    # 100000000000000000000000000319

i 番目の素数　prime(::Type{<:Integer}=Int, i::Integer)

100000 番目の素数は 1299709

prime(100000) # 1299709

primes() で生成される 1299709 までの素数の個数を数えれば，100000 になっていることが確認できる。

length(primes(1299709)) # 100000

素数判定　isprime(n::Integer) -> Bool

isprime(100000000000000000000000000319) # true

BigInt の場合は，確率的素数判定を行う。偽陽性率（素数ではないのに素数であると判定する率）は 0.25^reps。デフォルトでは reps = 25 なので，0.25^25 = 8.881784197001252e-16 に設定される。
isprime(x::BigInt, [reps = 25]) -> Bool

isprime(big(100000000000000000000000000319)) # true
isprime(big"100000000000000000000000000319") # true

isprime(big(100000000000000000000000000318)) # false

メルセンヌ素数判定　ismersenneprime(M::Integer; [check::Bool = true]) -> Bool

メルセンヌ数とは、2の冪よりも 1 小さい自然数、すなわち 2^n − 1（n は自然数）の形の自然数のことである。
これが素数の場合はメルセンヌ素数と呼ばれる。

ismersenneprime(2^11 - 1) # false
ismersenneprime(2^13 - 1) # true

素数の篩　primesmask([lo,], hi)

lo から hi までの数に対して，素数なら 1，合成数なら 0 の BitVector を返す。

mask = primesmask(10)
#=
10-element BitVector:
 0
 1
 1
 0
 1
 0
 1
 0
 0
 0
=#

mask が 1 である数を取り出す。

(1:10)[mask]
#=
4-element Vector{Int64}:
 2
 3
 5
 7
=#

ラディカル　radical(n::Integer)

偶数個の素因数を持たないように選択した素因数の累積。

n = 2*2*3*7*7*7 # 4116
radical(n) # 42 = 2 * 3 * 7
factor(radical(n)) # 2 * 3 * 7

トーシエント数　totient(f::Factorization{T}) -> T
                   totient(n::Integer) -> Integer

オイラーのトーシエント関数（ ϕ(n); totient function) の値
m ≦ n で，gcd(m, n) == 1 つまり，m, n が互いに素であるものの数

n が　素数の場合は 1 〜 n-1 と n は互いに素であるので，ϕ(11) は 10

sum(gcd(m, 11) == 1 for m = 1:11)
totient(11) # 10

n が合成数の場合は，n - 1 より小さな値になる

sum(gcd(m, 100) == 1 for m = 1:100) # 40
totient(100) # 40

Julia の小ネタ--026 R の rle(), Run Length Encoding

同じ文字列が何個続くか勘定する。

aabbbcccc なら，a が 2個，b が 3個，c が 4個なので，([2, 3, 4], ["a", "b", "c"]) のようなタプルを返す。

function rle(x)
    """ Run Length Encoding
    Compute the lengths and values of runs of equal values in a vector – or the reverse operation.
    
    rle("001111000000011")
    # ([2, 4, 7, 2], ["0", "1", "0", "1"])

    rle("aabcccbbbaab")
    # ([2, 1, 3, 3, 2, 1], ["a", "b", "c", "b", "a", "b"])
    """
    s = split(x, "")
    n = length(x)
    n == 0 && return (Int[], "")
    y = s[2:end] .!= s[1:end-1]
    i = vcat((1:length(y))[y], n)
    diff(vcat(0, i)), string.(s[i])rle(
end

rle("001111000000011") # ([2, 4, 7, 2], ["0", "1", "0", "1"])
rle("aabcccbbbaab")    # ([2, 1, 3, 3, 2, 1], ["a", "b", "c", "b", "a", "b"])

StatsBase にも rle() という関数がある。全く同じことをやるのだが，引数はベクトルである。返すタプルも，(種類, 個数) の順。

 

Julia で多項式，求解，多項式近似

using Polynomials
using Plots
gr(xlabel="\$x\$", ylabel="\$y\$", size=(400, 300))

多項式の定義，多項式の値，多項式の解

Polynomial() の目的は，1つには与えた係数ベクトルから「多項式 Polynomial{Float64, :x}」を作ることである。
係数ベクトルは定数項（0 次）から，1，2，3 次... の順に指定する。該当次数がないときは 0 を指定する。

p = Polynomial([1, 0, 3, 4])　# 1 + 3∙x^2 + 4∙x^3
typeof(p)                     # Polynomial{Float64, :x}

戻り値 p は関数でもあるので，plot() の第1引数（関数）としても使える。

plot(p, tick_direction=:out, title="\$1 + 3x^2 + 4x^3\$", label="")
hline!([0])
savefig("fig1.png")

また，関数なので当然ながら，指定した値に対応する多項式の値を求めることもできる。

p(-1)    # 0
p(-0.75) # 1.0
p(0)     # 1

特に，多項式=0 の解を求めるには，roots() を使う。複素解も求まる。

roots(p)
#=
3-element Vector{ComplexF64}:
                -1.0 + 0.0im
 0.12500000000000006 - 0.4841229182759272im
 0.12500000000000006 + 0.4841229182759272im
=#

多項式は，多項式 = 0 の解を a1, a2, ...としたとき，解のベクトル [a1, a2, ...] を指定して fromroots() で作成することもできる。

(x + 1.5) * (x - 2.2) * (x - 3.7) = 0 のとき，以下のように与える。

p =  fromroots([-1.5, 2.2, 3.7]) # 12.21 - 0.7099999999999995∙x - 4.4∙x^2 + 1.0∙x^3
plot(p, tick_direction=:out, title="\$12.21 - 0.71^x - 4.4x^2 + 1.0 x^3\$", label="")
hline!([0])
savefig("fig2.png")

p(-1.5) # -3.552713678800501e-15
p(2.2)  # -1.7763568394002505e-15
p(3.7)  # 0.0

Machine.double.epsilon は eps() = 2.220446049250313e-16 なので，8eps()~16eps() 程度の誤差である。

eps()   # 2.220446049250313e-16
8eps()  # 1.7763568394002505e-15
16eps() # 3.552713678800501e-15

roots(p)
#=
3-element Vector{Float64}:
 -1.4999999999999998
  2.200000000000002
  3.7
=#
p(-1.4999999999999998) # 1.7763568394002505e-15
p(2.200000000000002)   # -1.2434497875801753e-14
p(3.7)                 # 0.0

多項式に用いる変数の指定

第二引数に「シンボル」で与えることができる（デフォルトは，前述の通り :x とするのと同じ）。

r = Polynomial([1,2,3], :a) # 1 + 2∙a + 3∙a^2

多項式の四則演算

四則演算の一部が使える。ただし，多項式に使われる変数が同じでなければならない。

p = Polynomial([1, 2])     # 1 + 2∙x
q = Polynomial([1, 0, -1]) # 1 - x^2

定数倍

2p  # 2 + 4∙x
q/2 # 0.5 - 0.5∙x^2

定数和

3 + p # 4 + 2∙x

和

p +4q # 5 + 2∙x - 4∙x^2

差

p - q # 2∙x + x^2

積

p * q # 1 + 2∙x - x^2 - 2∙x^3

割り算

q ÷ p        # 0.25 - 0.5∙x
div(q, p)    # 0.25 - 0.5∙x
rem(q, p)    # 0.75
divrem(q, p) # (Polynomial(0.25 - 0.5*x), Polynomial(0.75))

不定積分

y = Polynomial([1, 0, -1]) # 1 - x^2
z = integrate(y)           # 1.0∙x - 0.3333333333333333∙x^3

微分

derivative(z) # 1.0 - 1.0∙x2　--> y に戻る

y = (x + 1) * (x - 2) * (x -4) の極値を求める。

多項式の定義

y = fromroots([-1, 2, 4])
plot(y, tick_direction=:out, title="\$8 + 2x - 5x^2 + x^3\$", label="")
savefig("fig3.png")

y の導関数を求めて(derivative)，0 になるときの値 x を求める(roots)。

x = roots(derivative(y))
#=
2-element Vector{Float64}:
 0.2137003521531088
 3.1196329811802244
=#

vline!([x])
savefig("fig4.png")

多項式近似（多項式補間）

前報で述べたように，滑らかな単調増加・単調減少関数は多項式によって有効に近似できる。

正規分布の上側確率を求める近似式を作ってみる。

using Rmath

z値 に対する p 値の 9 対のデータに基づく。

z = 1:0.5:5          # z 値, 1.0:0.5:5.0
p = pnorm.(z, false) # p 値
#=
9-element Vector{Float64}:
 0.15865525393145705
 0.06680720126885807
 0.022750131948179212
 0.006209665325776135
 0.0013498980316300946
 0.00023262907903552494
 3.1671241833119924e-5
 3.39767312473006e-6
 2.866515718791939e-7
=#

理論的に 8 次の近似式まで可能。

scatter(z, p, tick_direction=:out, title="pnorm(x)", label="sample point")
approx = fit(z, p); # length(z) - 1 次式（8 次式）で近似

見た目では差がわからないほどよく近似されているように見えるが。

plot!(approx, [1,5]..., label="approximation")
savefig("fig5.png")

近似式を求めたのと同じ範囲の z をさらに細かく(0.01刻み）にして，真値と近似値の差をプロットしてみる。

z2 = collect(1:0.01:5);
truevalue = similar(z2);
approximatevalue = similar(z2);
for (i, z3) in enumerate(z2)
    truevalue[i] = pnorm(z3, false)
    approximatevalue[i] = approx(z3)
end

少なくとも小数点以下 4 桁程度まで正確であることがわかる。

plot(z2, (truevalue .- approximatevalue),
    xlabel="\$z\$", ylabel="truevalue - approximatevalue", label="")
savefig("fig6.png")

Julia で多項式回帰，Polynomials

x = range(0, 5, length=6)
y = @. exp(-x)
#=
6-element Vector{Float64}:
 1.0
 0.36787944117144233
 0.1353352832366127
 0.049787068367863944
 0.01831563888873418
 0.006737946999085467
 =#

というデータがある。y を x の多項式 a0 + a1*x + a2*x^2 + ... で予測すること（多項式回帰）を考える。

多項式回帰は x^2, x^3, ... を前もって用意しておけば重回帰プログラムを使えばできる。
理論上は x の次数は length(y) - 1 までなので，例の場合は 5 次の項まで使える。

using DataFrames
df = DataFrame(:x => x, :x2 => x.^2, :x3 => x.^3, :x4 => x.^4, :x5 => x.^5, :y => y)
#=
6 rows × 6 columns

 x x2 x3 x4 x5 y
 Float64 Float64 Float64 Float64 Float64 Float64
1 0.0 0.0 0.0 0.0 0.0 1.0
2 1.0 1.0 1.0 1.0 1.0 0.367879
3 2.0 4.0 8.0 16.0 32.0 0.135335
4 3.0 9.0 27.0 81.0 243.0 0.0497871
5 4.0 16.0 64.0 256.0 1024.0 0.0183156
6 5.0 25.0 125.0 625.0 3125.0 0.00673795
=#

重回帰プログラムは GLM にある lm() を使う。

using GLM

lmans2 = lm(@formula(y ~ x + x2), df); # 2 次式
coef(lmans2)
#= 0, 1, 2 次の係数（以下同様）
3-element Vector{Float64}:
  0.9313226293358078
 -0.5231411796668056
  0.0697679508663814
=#

lmans3 = lm(@formula(y ~ x + x2 + x3), df); # 3 次式
coef(lmans3)
#=
4-element Vector{Float64}:
  0.9917995957122201
 -0.7993193261190775
  0.22096036680740455
 -0.0201589887921363
=#

lmans4 = lm(@formula(y ~ x + x2 + x3 + x4), df); # 4 次式
coef(lmans4)
#=
5-element Vector{Float64}:
  0.9995995032132226
 -0.9293177844691326
  0.36070870953370876
 -0.0656584492146525
  0.0045499460422515035
=#

lmans5 = lm(@formula(y ~ x + x2 + x3 + x4 + x5), df); # 5 次式
coef(lmans5)
#=
6-element Vector{Float64}:
  0.9999999999998723
 -0.9762026083014735
  0.4413086878624456
 -0.11144858183069759
  0.015062986693944548
 -0.0008410432521384522
=#

5次多項式による予測値

predict(lmans5)
#=
6-element Vector{Float64}:
 0.9999999999998723
 0.36787944117195287
 0.13533528323580923
 0.04978706836849145
 0.018315638888491637
 0.006737946999125999
=#

多項式回帰を行う fit() が Polynomials にある。
次数が高い項を使うと，正規方程式は不安定になるが，Polynomials.fit() だと，不安定性を軽減する計算アルゴリズムが採用される。

using Plots, Polynomials

注：以下では Polynomials.fit() として使っているが，これは先に使った GLM の fit() と識別するためなので，通常は fit() だけでよい。

f2 = Polynomials.fit(x, y, 2); # 2次式で予測
# 0.9313226293358072 - 0.5231411796668045∙x + 0.06976795086638117∙x2

f3 = Polynomials.fit(x, y, 3); # 3次式で予測
# 0.9917995957122127 - 0.7993193261190567∙x + 0.22096036680739517∙x2 - 0.02015898879213521∙x3

f4 = Polynomials.fit(x, y, 4); # 4次式で予測
# 0.9995995032131947 - 0.9293177844687606∙x + 0.360708709533328∙x2 - 0.06565844921453219∙x3 + 0.004549946042239714∙x4

f5 = Polynomials.fit(x, y); # 最高次数 = length(x) - 1 = 5次式で予測
# 1.0 - 0.9762026083107394∙x + 0.4413086878778398∙x2 - 0.11144858183923873∙x3 + 0.015062986695871163∙x4 - 0.0008410432522905108∙x5

データにある x の範囲内だと次数が高いほど予測が上手くできるように見えるが，

scatter(x, y, grid=false, tick_direction=:out,
        markerstrokewidth=0, size=(400, 300),
        xlabel="\$x\$", ylabel="\$y\$", label="Data")
plot!(f2, extrema(x)..., label="Order2")
plot!(f3, extrema(x)..., label="Order3")
plot!(f4, extrema(x)..., label="Order4")
plot!(f5, extrema(x)..., label="Order5")
savefig("fig1.png")

範囲外では全くダメなのがよく分かる。

scatter(x, y, grid=false, tick_direction=:out,
        markerstrokewidth=0, size=(400, 300),
        xlims=(-3,8),
        xlabel="\$x\$", ylabel="\$y\$", label="Data")
plot!(f2, [-3,8]..., label="Order2")
plot!(f3, [-3,8]..., label="Order3")
plot!(f4, [-3,8]..., label="Order4")
plot!(f5, [-3,8]..., label="Order5")
savefig("fig2.png")

また，単調減少や単調増加の場合は目立たないが，そうでない場合にはデータ範囲内でもデータ点は通っても，それ以外の所は凸凹で，予測とはとてもいえないことが明らかである。

x2 = 0:5
y2 = [2,6,9,4,11,7]
f52 = Polynomials.fit(x2, y2); # degree = length(x) - 1 = 5
scatter(x2, y2, grid=false, tick_direction=:out,
        markerstrokewidth=0, size=(400, 300),
        xlabel="\$x\$", ylabel="\$y\$", label="")
plot!(f52, [-0.3,5.3]..., label="")
savefig("fig3.png")

Julia に翻訳--240 陽性反応適中率(ppv)の差の検定，陰性反応適中率(npv)の差の検定

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

陽性反応適中率の差の検定
http://aoki2.si.gunma-u.ac.jp/JavaScript/ppv.html

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

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

==========#

using Rmath, Printf

function ppv(a, b, c, d, e, f, g, h)
  ppv(1, a, b, c, d, e, f, g, h)
  ppv(2, d, c, b, a, h, g, f, e)
end

function ppv(sw, a, b, c, d, e, f, g, h)
  name = ["ppv", "npv"]
  title = ["陽性反応適中率", "陰性反応適中率"]
  println("*****  $(title[sw])($(name[sw]))の差の検定  *****\n")
  if (sw == 1)
    @printf("%s1 = %g\n", name[sw], (a+b)/(a+b+e+f))
    @printf("%s2 = %g\n", name[sw], (a+c)/(a+c+e+g))
  else
    @printf("%s1 = %g\n", name[sw], (e+f)/(a+b+e+f))
    @printf("%s2 = %g\n", name[sw], (e+g)/(a+c+e+g))
  end
  mi = 2(a+e)+b+c+f+g
  zbar = (a+c+e+g)/mi
  dbar = (2a+b+c)/mi
  stat = (a * (1 - 2zbar) - b * zbar + c * (1 - zbar))^2 /
    ((1-dbar)^2 * (a * (1 - 2zbar)^2 + b*zbar^2 + c * (1 - zbar)^2) + dbar^2 *
     (e * (1 - 2zbar)^2 + f * zbar^2 + g * (1 - zbar)^2))
  p = pchisq(stat, 1, false)
  @printf("T%s = %g\n", name[sw], stat)
  @printf("P value = %.5f (%g)\n\n", p, p)
end

ppv(786, 183, 29, 25, 69, 176, 46, 151)

#==
*****  陽性反応適中率(ppv)の差の検定  *****

ppv1 = 0.798188
ppv2 = 0.876344
Tppv = 47.5813
P value = 0.00000 (5.2769e-12)

*****  陰性反応適中率(npv)の差の検定  *****

npv1 = 0.784861
npv2 = 0.611215
Tnpv = 39.5222
P value = 0.00000 (3.24349e-10)
==#