ノートブック Convert ipynb to html, tex, pdf
を実行すると, このファイルから jmd, html, tex, pdf ファイルが作成される. jmdファイルはJulia言語のコードを含むmarkdownファイルである.
$ \newcommand\real{\operatorname{Re}} \newcommand\imag{\operatorname{Im}} $
using Plots
gr()
##ENV["PLOTS_TEST"] = "true"
using SpecialFunctions
using Distributions
using SymPy
ガンマ函数とベータ函数は次のように定義される: $\real s, \real p, \real q > 0$ のとき
$$ \begin{aligned} & \Gamma(s) = \int_0^\infty e^{-x} x^{s-1}\,dx, \\ & B(p,q) = \int_0^1 x^{p-1} (1-x)^{q-1}\,dx \end{aligned} $$と定義される.
注意: 複数行の数式は二重のドルマークで囲んだ aligned モードを使うことにした. weave()
函数で作成した tex ファイルを訂正を施す必要がある.
ガンマ函数は階乗の連続的補間になっており, 急激に増大する函数になる. だから, グラフを描き易いように対数を取ったガンマ函数のグラフを描いてみよう. せっかくなので, グラフの中でStirlingの近似公式
$$ \log\Gamma(s) \approx s \log s - s - \frac{1}{2}\log s + \log\sqrt{2\pi} \quad\text{as $s\to\infty$} $$の右辺と比較してみよう. 以下のプロットを見ればわかるようにほぼぴったり一致している.
lstirling(s) = s*log(s) - s - log(s)/2 + log(√(2π))
s = range(0.1, 10, length=400)
plot(size=(500, 300), legend=:topleft)
plot!(s, lgamma.(s), label="log Gamma(s)", lw=2)
plot!(s, lstirling.(s), label="Stirling approximation", ls=:dash, lw=2)
$\theta,k,a,b>0$ に対して, ガンマ分布とベータ分布の確率密度函数がそれぞれ次のように定義される:
$$ \begin{aligned} & p(x) = \frac{1}{\theta^k \Gamma(k)} e^{-x/\theta} x^{k-1} \quad (x>0), \\ & q(x) = \frac{1}{B(a,b)} x^{a-1}(1-x)^{b-1} \quad (0 < x < 1). \end{aligned} $$これらは $k,q,b$ が大きなとき, それぞれ正規分布の確率密度函数
$$ \begin{aligned} f(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-(x-\mu)^2/(2\sigma^2)} \end{aligned} $$の $(\mu,\sigma^2)=(k\theta, k\theta^2), (a/(a+b), ab/((a+b)^2(a+b+1)))$ の場合でよく近似されるようになる(中心極限定理の特別な場合). そのことをグラフを描いて確認しよう.
k, θ = 100, 1/80
gdist = Gamma(k, θ)
ndist = Normal(k*θ, √k*θ)
x = range(0, 2.7, length=400)
plot(size=(500, 300))
title!("Normal approximation of Gamma distribution", titlefontsize=12)
plot!(legend=:topright)
plot!(x, pdf.(gdist, x), label="Gamma($k, $θ)", lw=2)
plot!(x, pdf.(ndist, x), label="Normal approximation", ls=:dash, lw=2)
a, b = 40, 60
bdist = Beta(a, b)
ndist = Normal(a/(a+b), √(a*b/((a+b)^2*(a+b+1))))
x = range(0, 1, length=400)
plot(size=(500, 300))
title!("Normal approximation of Beta distribution", titlefontsize=12)
plot!(legend=:topright)
plot!(x, pdf.(bdist, x), label="Beta($a, $b)", lw=2)
plot!(x, pdf.(ndist, x), label="Normal approximation", ls=:dash, lw=2)
上で使ったガンマ分布の確率密度函数 $p(x)=e^{-x/\theta}x^{k-1}/(\theta^k\Gamma(k))$ ($x>0$) が確率密度函数であることを示すためには, それの $0$ から $\infty$ までの積分が $1$ になること, すなわち,
$$ \int_0^\infty e^{-x/\theta} x^{k-1}\,dx = \theta^k \Gamma(k) \quad (\theta, k > 0) $$が成立することを示さなければいけない. この公式は左辺を $x=\theta y$ で置換することによって得られる.
そのことはSymPyを使っても確かめられる. (注意: 以下の計算はθをSymPyの変数として使用しているので Python 2.7 では不可能. Python 3.x を使いましょう.)
k, θ, x = symbols("k θ x", positive=true)
sol = simplify(integrate(exp(-x/θ)*x^(k-1), (x, 0, oo)))
solstr = replace(sympy.latex(sol), "θ"=>"\\theta")
display("text/html", raw"$$\int_0^\infty e^{-x/\theta} x^{k-1}\,dx = " * solstr * raw"$$");
注意: LaTeXStrings.jl の latexstring()
を使用して同じ数式を表示させると, weave(~, doctype="md2html")
でうまく数式が表示されなくなる. 上のように mimetype text/html
で display
を使用し, 行末に ;
を付ければうまく表示されるようになる.
Riemannのゼータ函数は $\real s > 1$ において
$$ \zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s} \quad (\real s > 1) $$と定義され, 唯一の極 $s=1$ を除いた複素平面全体に解析接続される.
Riemannのゼータ函数の $0\leqq \real s\leqq 1$ (critical strip)における零点を非自明な零点と呼ぶ.
Riemann予想: Riemannのゼータ函数の非自明な零点はすべて直線 $\real s=1/2$ 上に乗っている.
Riemannのゼータ函数の非自明な零点は素数の分布の精密な評価と関係している.
Riemann予想の成立を(数学的な証明にはならないが)数値的な計算で確認してみよう
f(s) = max(min(log(abs(zeta(s))), 3), -5)
x = range(0, 1, length=100)
y = range(10, 45, length=400)
s = @. x + im*y'
plot(size=(750, 180))
title!("log(abs(zeta(x+iy)))", titlefontsize=12)
heatmap!(y, x, f.(s), color=:thermal)
xlabel!("y")
ylabel!("x")
hline!([0.5], ls=:dot, color=:cyan, label="")
f(s) = abs(zeta(s))
x = 1/2
y = range(10, 45, length=1000)
s = @. x + im*y
plot(size=(720, 200))
title!("abs(zeta(0.5+iy)))", titlefontsize=12)
plot!(y, f.(s), label="")
xlabel!("y")
直線 $\real s = 1/2$ での $\zeta(s)$ の値は何度も繰り返し複素平面の原点 $0$ を通る.
x = 1/2
y = range(0, 45, length=1000)
s = @. x + im*y
w = zeta.(s)
plot(size=(500,400), legend=:topleft)
title!("zeta(0.5+it) for 0 < t < 45", titlefontsize=12)
plot!(real(w), imag(w), label="zeta(0.5+it)")
scatter!([0], [0], markersize=2, label="")
xlabel!("x")
ylabel!("y")
annotate!([(-1.25, 0.1, "zeta(0.5)", 9)])
annotate!([(2.55, 1.7, "zeta(0.5+45i)", 9)])
直線 $\real s = 0.6$ での $\zeta(s)$ の値を計算すると複素平面の原点を避けて通っていることがわかる.
x = 0.6
y = range(0, 45, length=1000)
s = @. x + im*y
w = zeta.(s)
plot(size=(500,400), legend=:topleft)
title!("zeta(0.6+it) for 0 < t < 45", titlefontsize=12)
plot!(real(w), imag(w), label="zeta(0.6+it)")
scatter!([0], [0], markersize=2, label="")
xlabel!("x")
ylabel!("y")
annotate!([(-1.6, 0.1, "zeta(0.6)", 9)])
annotate!([(2.55, 1.4, "zeta(0.6+45i)", 9)])