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