JuliaとJupyterのすすめ

注意: 筆者は「ライトユーザー」なので深い話はできない.

このファイルは以下の研究会のために準備された:

この文書は以下の場所で閲覧できる:

私によるJulia+Jupyterの環境構築の最近の記録については次のメモを参照せよ:

In [1]:
print("using Plots ->")
@time using Plots
gr(legend=false); ENV["PLOTS_TEST"] = "true"

print("\nplot(sin) ->")
@time plot(sin, size=(300, 160)) |> display # コンパイル

print("\nusing PyPlot ->")
@time using PyPlot: PyPlot, plt

print("\nusing DifferentialEquations ->")
@time using DifferentialEquations

using Base64
using Combinatorics
using Distributed
using Distributions
using Libdl
using LinearAlgebra
using Optim
using Primes
using ProgressMeter
using Random: seed!
using RCall
using SpecialFunctions
using Statistics
using SymPy: SymPy, sympy, Sym, @vars, @syms, simplify, oo, PI

ldisp(x...) = display("text/html", raw"$$" * prod(x) * raw"$$")

showimg(mime, fn) = open(fn) do f
    base64 = base64encode(f)
    display("text/html", """<img src="data:$mime;base64,$base64">""")
end
using Plots ->  6.497952 seconds (10.35 M allocations: 547.218 MiB, 5.71% gc time)
-4 -2 0 2 4 -1.0 -0.5 0.0 0.5 1.0
plot(sin) -> 25.666761 seconds (56.63 M allocations: 2.761 GiB, 6.50% gc time)

using PyPlot ->  7.997156 seconds (14.13 M allocations: 696.137 MiB, 5.28% gc time)

using DifferentialEquations -> 13.710626 seconds (35.83 M allocations: 2.302 GiB, 8.87% gc time)

R version 3.5.3 (2019-03-11) -- "Great Truth"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

Out[1]:
showimg (generic function with 1 method)

パッケージの読み込みと最初の実行時のコンパイルにそれなりに時間が採られることに注意.

Jupyter notebook では最初のセルで最初のコンパイルが実行されるようにしておき, そのセルを実行してから, プログラムのコードの入力を始めると時間の節約になる.

Jupyterノートブック

  • このファイルはJupyterノートブックの実例になっている.

  • JupyterではJulia, Python, R, Rubyなどのプログラミング言語のノートブックを作れる.

  • ノートブックにはブラウザ経由でアクセスする.

  • ノートブックには以下をまとめることができる:

    • プログラムのコード
    • その実行結果(グラフのプロットを含む)
    • 整形された数式を含む説明(markdown)
  • Jupterノートブックの表示にGitHubなどが対応している(nbviewerが非常に便利).

  • 試行錯誤に向いている.

  • 数百行程度のプログラミングにはこれで十分.

  • NbextensionsGist-itを使えば, ワンタッチでノートブックをGitHub Gistに保存してバージョン管理できる.

  • NbextensionsのRISEを使えば, Jupyterノートブックでプレゼン可能.

  • nbviewerのView as Slidesボタンもプレゼン可能(Slides).

JupyterはPythonのライブラリの1つ.

Anacondaを入れればJupyterも使えるようになっている.

追記(2019-05-26): Free Wolfram Engine が公開された. それとJupyterを組み合わせると「無料でMathematicaを使用」のような状態になる(実際にはJupyter側の対応が未成熟なので制限がある). 詳しくは

を参照して欲しい.

Python 2.7だと小2開発者に馬鹿にされるので注意

注意: この節はジョークのつもりで入れた.

整形された数式を含む説明の例

パラメーター $w$ 付きの $y$ に関する確率密度函数 $p(y|w)$ と事前分布と呼ばれる $w$ に関する確率密度函数 $\varphi(w)$ とサンプル $Y_1,\ldots,Y_n$ から得られるベイズ推測法における予測分布 $p^*(y)$ は以下のように定義される:

$$ p^*(y) = \frac{Z(Y_1,\ldots,Y_n,y|\varphi)}{Z(Y_1,\ldots,Y_n|\varphi)} = Z(y|\psi). $$

ここで

$$ \begin{aligned} & Z(y_1,\ldots,y_n|\varphi) = \int \prod_{k=1}^n p(Y_k|w)\cdot\varphi(w)\,dw, \\ & \psi(w) = \frac{\prod_{k=1}^n p(Y_k|w)\cdot\varphi(w)}{\int \prod_{k=1}^n p(Y_k|w)\cdot\varphi(w)\,dw}. \end{aligned} $$

$Z(y_1,\ldots,y_n|\varphi)$ は分配函数と呼ばれ, $Z(Y_1,\ldots,Y_n|\varphi)$ は周辺尤度と呼ばれ, $\psi(w)$ は事後分布と呼ばれる.

Julia言語

Julia言語は最新のプログラミング言語.

  • v1.0.0になったのは2018年8月. (2019-03-15現在はv1.1.0)

  • 無料. オープンソース.

  • JuliaBoxを使えばブラウザさえあれば使える.

  • スクリプト言語のように気軽に使える.

  • 計算が速い.

  • 私の場合には(個人差あり!), Cで書いてgccでコンパイルするよりも速いことが多い.

  • GNU OctavescilabのようなMATLABクローンのユーザーにはJulia言語に引っ越しすることを強く勧める.

  • Python, R, Rubyなどとは共存して行くことになるだろう.

私によるWindows環境でのJulia言語環境の作り方の記録:

私によるJulia言語の使用例のギャラリー:

v1.0がリリースされてから急激に人気が高まる

Heaviside函数の成分が追加されたのは2018年8月頃すなわちv1.0がリリースされたときである.

Just-in-time でコンパイル

Julia言語では

  • 函数を実行したときに,

  • 与えられた引数の型情報を使って,

  • その函数をネイティブコードにコンパイルしてから実行する.

こういう仕組みなので高速に計算したいコードは函数の中に入れなければいけない.

最初の実行時にはコンパイルにも時間が取られることに注意.

以下は円周率のモンテカルロ計算で比較.

In [2]:
### 函数にせずに実行

@time begin
    L = 10^7
    c = 0
    for i in 1:L
        global c
        c += ifelse(rand()^2+rand()^2  1, 1, 0)
    end
    4c/L
end
  1.418559 seconds (40.00 M allocations: 762.914 MiB, 10.46% gc time)
Out[2]:
3.1410228
In [3]:
### 函数にして実行

function simpi(L=10^7)
    c = 0
    for i in 1:L
        c += ifelse(rand()^2+rand()^2  1, 1, 0)
    end
    4c/L
end
simpi(10^5) # simpi(::Int64)がコンパイルされる

@time simpi()
  0.047364 seconds (36.00 k allocations: 1.967 MiB)
Out[3]:
3.141308

実行コードを函数の中に入れるだけで数十倍計算が速くなった.

Julia言語では各函数が引数の型情報に基いて just-in-time でコンパイルされる.

gccとの比較

円周率のモンテカルロ計算(1億回)で比較.

In [4]:
versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)
Environment:
  JULIA_CMDSTAN_HOME = C:\CmdStan
  JULIA_NUM_THREADS = 4
  JULIA_PKGDIR = C:\JuliaPkg
In [5]:
run(`gcc --version`); # ほとんど使っていないのでバージョンが古い
gcc (Rev3, Built by MSYS2 project) 6.3.0
Copyright (C) 2016 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

In [6]:
### gcc

C_code= raw"""
#include <time.h>
#include <stdlib.h>

double simpi(unsigned long n){
    double x,y;
    srand((unsigned)time(NULL));
    double R = (double) RAND_MAX;
    unsigned long count = 0;
    for(unsigned long j=0; j < n; j++){
        x = ((double) rand())/R;
        y = ((double) rand())/R;
        if(x*x + y*y <= 1){
            count++;
        }
    }
    return ((double) 4.0)*((double) count)/((double) n);
}
"""

filename = tempname()
filenamedl = filename * "." * Libdl.dlext

open(`gcc -Wall -O3 -march=native -xc -shared -o $filenamedl -`, "w") do f
     print(f, C_code)
end

## run(`ls -l $filenamedl`);

const simpi_gcc_rand_lib = filename
simpi_gcc_rand(n::Int64) = ccall((:simpi, simpi_gcc_rand_lib), Float64, (Int64,), n)
simpi_gcc_rand(10^5)

@time simpi_gcc_rand(10^8)
  2.537506 seconds (1.94 k allocations: 116.293 KiB)
Out[6]:
3.14153524
In [7]:
### Julia

@time simpi(10^8)
  0.392231 seconds (6 allocations: 192 bytes)
Out[7]:
3.14148068

gccの側が圧倒的に遅い. Juliaの方が計算が数倍速い.

これはgccのデフォルトの rand() が遅いから.

そもそもgccの rand() は擬似乱数の質の点でモンテカルロシミュレーションには向かない.

Julia言語は dSFMT(Double precision SIMD-oriented Fast Mersenne Twister) をデフォルトの擬似乱数発生器として使っている.

gccの側でもdSFMTを使えば, Julia言語より少し速くなる.

注意: MT19937 より dSFMT の方が擬似乱数の質も上でかつ擬似乱数の発生速度も約3倍程速い(検証用コード). この点に MT19937 (所謂メルセンヌツイスター)の使用者は注意した方がよい. ライブラリ選択に関するこの手の問題は無数にあるので, 私程度の知識と実力だと gcc を使った場合よりも Julia を使った場合の方が計算が速くなることが多い. Julia言語は信頼できてかつ高速なライブラリの寄せ集めにもなっている.

以下で円周率のモンテカルロ計算(10億回)で比較.

In [8]:
### gcc with dSFMT

## Download http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/dSFMT-src-2.2.3.tar.gz
##
## Extract
##   dSFMT-src-2.2.3/dSFMT.h
##   dSFMT-src-2.2.3/dSFMT.c

C_code= raw"""
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define HAVE_SSE2
#define DSFMT_MEXP 521
#include "dSFMT-src-2.2.3/dSFMT.h"
#include "dSFMT-src-2.2.3/dSFMT.c"

double simpi(unsigned long n){
    srand((unsigned)time(NULL));
    dsfmt_t dsfmt;
    dsfmt_init_gen_rand(&dsfmt, rand());
    unsigned long count = 0;
    double x,y;
    for(unsigned long j = 0; j < n; j++){
        x = dsfmt_genrand_close_open(&dsfmt);
        y = dsfmt_genrand_close_open(&dsfmt);
        if(x*x + y*y <= 1){
            count++;
        }
    }
    return ((double)4.0)*((double)count)/((double)n);
}
"""

filename = tempname()
filenamedl = filename * "." * Libdl.dlext

open(`gcc -Wall -O3 -march=native -xc -shared -o $filenamedl -`, "w") do f
     print(f, C_code)
end

## run(`ls -l $filenamedl`);

const simpi_gcc_dSFMT_lib = filename
simpi_gcc_dSFMT(n::Int64) = ccall((:simpi, simpi_gcc_dSFMT_lib), Float64, (Int64,), n)
simpi_gcc_dSFMT(10^5)

@time simpi_gcc_dSFMT(10^9)
  2.928520 seconds (1.41 k allocations: 84.691 KiB)
Out[8]:
3.14162012
In [9]:
### Julia

@time simpi(10^9)
  3.835884 seconds (6 allocations: 192 bytes)
Out[9]:
3.141565032

まとめ

  • gccを使った方が確かに計算は速くなったが, そのためには適切なライブラリの選定とダウンロードと使用が必要になった.

  • 上の例ではその作業は簡単だったが, もっと実戦的な例では非常に面倒なことになる.

  • Julia言語では多数の優れたライブラリがデフォルトで使えるようになっている.

  • Julia言語を使うとストレスが大幅に減る!

グルー言語としてのJulia

Julia言語は複数のプログラミング言語の「貼り合わせ」に使えるグルー言語(糊言語)でもある.

以上を見ればわかるように

  • Juliaではgccでコンパイルされた函数を簡単に使える.

この他にも

  • Fortranの函数も使える.

  • Pythonとの連携は極めて密(PyCall.jl).

  • Rとも連携(RCall.jl).

  • ……(以下略)

JuliaとPythonとRとのあいだのフレンドリーな連携がすでに行われているので, 今後も Julia, Python, R の3つは互いに影響を与えながら発展して行くものと期待される.

注意(2019-05-04更新): C++との連携も可能である(Cxx.jl). ただし, 2019-05-04の時点でWindows対応は初期段階に過ぎず, 使える機能と使えない機能がある. Windows環境であっても, g++でコンパイルして作ったshared libraryをJuliaで利用することができ, Julia内から直接C++を使用することもある 程度可能である. 詳しくは次のリンク先を参照:

Pythonのmatplotlibの利用

一般にグラフのプロットの仕方の学習コストは高い.

Pythonに慣れている人はmatplotlibをほぼそのままJuliaでも利用できる.

Julia言語におけるプトッロのためのライブラリとして次の2つがおすすめ:

私は片方に統一したりせずに両方を使っている. Plot.jlpyplot() バックエンドで使うこともよくある. それら以外にも

も安定して使用可能にできる環境ならば便利なようである.

Riemannのゼータ函数の絶対値のクリティカルストリップでのプロット

In [10]:
f(s) = max(min(log(abs(zeta(s))), 2), -4)
x = range(0, 1, length=100)
y = range(10, 45, length=400)
s = @. x + im*y'
@time w = f.(s)
plt.figure(figsize=(8,1.5))
plt.pcolormesh(imag(s), real(s), w, cmap="gist_rainbow_r")
plt.hlines(0.5, minimum(y), maximum(y), colors="grey", lw=0.5, linestyles="dashed")
plt.xlabel("y")
plt.ylabel("x")
plt.title("log(abs(zeta(x+iy)))");
  0.692026 seconds (1.49 M allocations: 75.519 MiB, 5.28% gc time)

Pythonのsympyの利用

Pythonで書かれたsympyという数式処理系をJuliaでも利用できる.

私によるJuliaにおけるSymPyの利用例

In [11]:
@vars a r positive=true
@vars x real=true
I = sympy.Integral(exp(-a*x^r), (x,0,oo))
sol = simplify(I.doit())
ldisp(sympy.latex(I), " = ", sympy.latex(sol))
$$\int_{0}^{\infty} e^{- a x^{r}}\, dx = a^{- \frac{1}{r}} \Gamma\left(1 + \frac{1}{r}\right)$$

数値計算用のコードの出力をSymPyを使って数式で表示できる!

In [12]:
### これは x² = ax + b を満たす x の連分数展開

function f(a, b, x, n)
    s = a + b/x
    for i in n:-1:1
        s = a+b/s
    end
    s
end

### 連分数で x² = x + 1 の解 x (黄金分割比)を数値計算

f(1,1,1,37), (1+√5)/2
Out[12]:
(1.618033988749895, 1.618033988749895)
In [13]:
### SymPyを使えば連分数を整形して表示できる.

@vars a b x
for n in 0:3
  ldisp("f(a,b,x,$n) = ", sympy.latex(f(a, b, x, n)))
end
$$f(a,b,x,0) = a + \frac{b}{x}$$
$$f(a,b,x,1) = a + \frac{b}{a + \frac{b}{x}}$$
$$f(a,b,x,2) = a + \frac{b}{a + \frac{b}{a + \frac{b}{x}}}$$
$$f(a,b,x,3) = a + \frac{b}{a + \frac{b}{a + \frac{b}{a + \frac{b}{x}}}}$$
In [14]:
R"""
library(ggplot2)
str(iris)
""";

R"""
p <- ggplot(iris, aes(x=Sepal.Length, y=Sepal.Width))
p + geom_point(colour="gray50", size=3) + geom_point(aes(colour=Species))
"""
'data.frame':	150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
Out[14]:
RObject{VecSxp}

Julia側のデータをR側で分析

In [15]:
n = 100
X = randn(n,2)
b = [2.0, 3.0]
y = X * b + randn(n,1)

@rput y
@rput X

R"mod <- lm(y ~ X-1)"
R"summary(mod)"
Out[15]:
RObject{VecSxp}

Call:
lm(formula = y ~ X - 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5326 -0.5529 -0.1147  0.3792  2.3786 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
X1  2.08503    0.08504   24.52   <2e-16 ***
X2  3.09001    0.08537   36.20   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9142 on 98 degrees of freedom
Multiple R-squared:  0.9473,	Adjusted R-squared:  0.9463 
F-statistic: 881.6 on 2 and 98 DF,  p-value: < 2.2e-16

ギャラリー

Plots.jlによるアニメーション

以下のようにして簡単かつそれなりに短時間でアニメーションも作成できる.

In [16]:
# 線分を描く函数
segment(A, B; color="black", kwargs...) = plot([A[1], B[1]], [A[2], B[2]]; color=color, kwargs...)
segment!(A, B; color="black", kwargs...) = plot!([A[1], B[1]], [A[2], B[2]]; color=color, kwargs...)
segment!(p, A, B; color="black", kwargs...) = plot!(p, [A[1], B[1]], [A[2], B[2]]; color=color, kwargs...)

# 円周を描く函数
function circle(O, r; a=0, b=2π, color="black", kwargs...)
    t = linspace(a, b, 1001)
    x(t) = O[1] + r*cos(t)
    y(t) = O[1] + r*sin(t)
    plot(x.(t), y.(t); color=color, kwargs...)
end
function circle!(O, r; a=0, b=2π, color="black", kwargs...)
    t = linspace(a, b, 1001)
    x(t) = O[1] + r*cos(t)
    y(t) = O[2] + r*sin(t)
    plot!(x.(t), y.(t); color=color, kwargs...)
end
function circle!(p, O, r; a=0, b=2π, color="black", kwargs...)
    t = linspace(a, b, 1001)
    x(t) = O[1] + r*cos(t)
    y(t) = O[1] + r*sin(t)
    plot!(p, x.(t), y.(t); color=color, kwargs...)
end
Out[16]:
circle! (generic function with 2 methods)
In [17]:
function MLIanim()
    lw = 1.0 # 太さ
    A, B, C, D = [1.5, 0], [3.5, 0], [5.5, 0], [7.5, 0]
    N = 9
    θ₀ = 2π/(2N)
    R = [
        cos(θ₀) -sin(θ₀)
        sin(θ₀)  cos(θ₀)
    ]
    V(θ) = [cos(θ), sin(θ)]
    r = 1.55*2π/(2N)

    aa = [0.5:0.025:1.5; 1.475:-0.025:0.525]
    prog = Progress(length(aa),1)
    anim = @animate for a in aa
        θ = a*2π/4
        p = plot(xlim=(-9, 9), ylim=(-9, 9))
        plot!(grid=false, legend=false, xaxis=false, yaxis=false)
        for k in 1:10
            RR = R^(2k-1)
            segment!(RR*A, RR*B,               lw=lw, color=:black)
            segment!(RR*B, RR*C,               lw=lw, color=:blue)
            segment!(RR*C, RR*D,               lw=lw, color=:red)
            segment!(RR*A, RR*(A+r*V(θ)),      lw=lw, color=:black)
            segment!(RR*A, RR*(A+r*V(-θ)),     lw=lw, color=:black)
            segment!(RR*B, RR*(B+1.5r*V(π-θ)), lw=lw, color=:black)
            segment!(RR*B, RR*(B+1.5r*V(π+θ)), lw=lw, color=:black)
            segment!(RR*C, RR*(C+2.0r*V(θ)),   lw=lw, color=:black)
            segment!(RR*C, RR*(C+2.0r*V(-θ)),  lw=lw, color=:black)
            segment!(RR*D, RR*(D+2.5r*V(π-θ)), lw=lw, color=:black)
            segment!(RR*D, RR*(D+2.5r*V(π+θ)), lw=lw, color=:black)
        end
        plot(p, size=(500, 500))
        next!(prog)
    end
    anim
end
Out[17]:
MLIanim (generic function with 1 method)

アニメーション

In [18]:
anim = MLIanim()
gifname = "dynamic_Muller-Lyer.gif"
gif(anim, gifname, fps = 15)
sleep(0.1)
showimg("image/gif", gifname)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:28
┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\msfd28\dynamic_Muller-Lyer.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\47Tik\src\animation.jl:90

特異モデルに近い場合の尤度函数

確率モデル(パラメーター $a,b$ 付きの $y$ に関する確率分布):

$$ \begin{aligned} & p(y|a,b) = a N(y) + (1-a)N(y-b), \\ & N(y) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}. \end{aligned} $$

サンプルを生成する真のモデル: $p(y|a_0, b_0)$.

$a_0=1$ または $b_0=0$ のとき上の確率モデルは特異モデルになる.

以下では特異モデルになってしまう $b_0=0$ に近い $(a_0,b_0)=(0.5, 0.1)$ の場合の尤度函数をプロットしてみる.

In [19]:
mixnormal(a,b) = MixtureModel([Normal(0,1), Normal(b,1)], [a, 1-a])
lpdf(a, b, y) = log(a+(1-a)*exp(b*y-b^2/2)) - y^2/2 - log(2π)/2
loglik(a, b, Y) = sum(lpdf(a, b, Y[i]) for i in 1:lastindex(Y))

function plot_lik(a₀, b₀, n; seed = 4649, bmin=-1.0, bmax=1.0)
    seed!(seed)
    Y = rand(mixnormal(a₀, b₀), n)
    L = 201
    a = range(0, 1, length=L)
    b = range(bmin, bmax, length=L)
    f(a, b) = loglik(a, b, Y)
    z = f.(a', b)
    zmax, k = findmax(z)
    #i, j = (k - 1) ÷ L + 1, mod1(k, L) # for Julia v0.6
    j, i = k.I
    z .= exp.(z .- zmax)

    sleep(0.1)
    plt.figure(figsize=(5,5))
    plt.pcolormesh(a, b, z, cmap="CMRmap")
    plt.scatter([a₀], [b₀], label="true", color="cyan")
    plt.scatter([a[i]], [b[j]], label="MLE", color="red")
    plt.legend()
    plt.xlim(0,1)
    plt.xlabel("a")
    plt.ylabel("b")
    plt.title("\$a_0 = $a₀,  b_0 = $b₀,  n = $n\$")
end
Out[19]:
plot_lik (generic function with 1 method)

サンプルサイズ 512 の場合

尤度函数は特異点集合 $a=1$ または $b=0$ に沿って拡がった形になっている.

最尤法(MLE)によるパラメーターの推定値が真の値 $(a_0,b_0)=(0.5,0.1)$ とは全然違う値になっている.

In [20]:
@time plot_lik(0.5, 0.1, 2^9);
  2.259230 seconds (2.65 M allocations: 130.699 MiB, 3.25% gc time)

サンプルサイズ 4096 の場合

さらにサンプルサイズを8倍に増やしても, 尤度函数は特異点集合 $a=1$ または $b=0$ に沿って拡がった形のままである.

やはり, 最尤法(MLE)によるパラメーターの推定値が真の値 $(a_0,b_0)=(0.5,0.1)$ とは全然違う値になっている.

サンプルサイズを $2^{18}=262144$ まで増やしても, 尤度函数は特異点集合に沿って広がったままで, 最尤法の結果も真の値に近付かないことを数値的に確認できる.

In [21]:
@time plot_lik(0.5, 0.1, 2^12);
  8.352434 seconds (96.16 k allocations: 3.204 MiB)

数行でランダムウォークをプロットできる

Bernoulli試行のShannon情報量の推定

In [22]:
### Shannon情報量の収束の様子
seed!(2019)
ber = Bernoulli(1/3)
S = entropy(ber)
nmax, ntrials = 2^10, 5
n = 1:nmax
a = cumsum(rand(ber, nmax, ntrials), dims=1)./n
y = @. -(a*log(ber.p) + (1-a)*log(1-ber.p))
plot(size=(500, 300), y)
hline!([S], color=:black, ls=:dash)
Out[22]:
0 250 500 750 1000 0.4 0.6 0.8 1.0