Juliaで作って学ぶベイズ統計学のJulia コードを実行してみた結果を，GItHubで公開してみています．
Jupyter Notebookが完成し次第，随時アップロードしていく予定です．
2022年9月14日水曜日
2022年9月10日土曜日
JuliaをアップデートしてJupiter Notebookでも使えるようにする
アップデートなので，すでにJuliaがインストールされている前提です．
Homebrewを使います．まずはターミナルから，以下のようしてJuliaをアップデートします．
% brew install julia --cask
以前は，以下のように cask install juliaだったのですが，エラーが出てしまいました．
% brew cask install julia
Error: `brew cask` is no longer a `brew` command. Use `brew <command> --cask` instead.
% brew cask install julia
アップデートが終わったら，Applicationフォルダに新しいJulia（以下の例では，最初に1.0がインストールされていて，新たに1.8をインストールした）のアイコンがあるので，アイコンをダブルクリックして立ち上げます．
その上で，"]"キーを押してパッケージモードに切り替えて，以下のように（add IJulia）IJuliaパッケージを追加します．
/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia ; exit;
% /Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia ; exit;
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.8.1 (2022-09-06)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
(@v1.8) pkg> add IJulia
この↑例では省略していますが，add IJuliaを実行すると，ずらずらとインストールが始まります．
IJuliaパッケージの追加が終わった段階で，Available kernelsを確認すると，julia-1.8 が現れます．
% jupyter kernelspec list
Available kernels:
julia-1.0 /Users/hide/Library/Jupyter/kernels/julia-1.0
julia-1.8 /Users/hide/Library/Jupyter/kernels/julia-1.8
python3 /usr/local/share/jupyter/kernels/python3
Jupyter Notebookを立ち上げて，インストールしたJuliaを使えるかどうか確認しておきます．
Julia（今回は1.8.1）が使えるようになっていれば成功です．
使わないカーネルや，存在しないカーネル（今回の例では，Julia 1.0.5）を削除するにはターミナルから以下のように入力します．なお，"KERNEL_NAME"には使わなくなった，もしくは存在しないカーネルの名前を書きます．
% jupyter kernelspec uninstall KERNEL_NAME
今回の例で，julia-1.0を削除すると以下のようになります．
% jupyter kernelspec uninstall julia-1.0
Kernel specs to remove:
julia-1.0 /Users/xxx/Library/Jupyter/kernels/julia-1.0
Remove 1 kernel specs [y/N]: y
[RemoveKernelSpec] Removed /Users/xxx/Library/Jupyter/kernels/julia-1.0
今回の例では，Juliaを1.0から1.8にアップデートしました．1.0でインストールしたパッケージは継承されないようなので，改めて追加する必要があります．以下の例は，とりあえず追加しておいたパッケージ達です．
(@v1.8) pkg> add CSV
(@v1.8) pkg> add Combinatorics
(@v1.8) pkg> add DataFrames
(@v1.8) pkg> add Distributions
(@v1.8) pkg> add JuMP
(@v1.8) pkg> add Plots
(@v1.8) pkg> add PyCall
(@v1.8) pkg> add PyPlot
(@v1.8) pkg> add RDatasets
(@v1.8) pkg> add StatsBase
(@v1.8) pkg> add StatsPlots
(@v1.8) pkg> add SumOfSquares
(@v1.8) pkg> add Statistics
(@v1.8) pkg> add Gadfly
(@v1.8) pkg> add LinearAlgebra
(@v1.8) pkg> add Optim
インストールしたパッケージを確認するには，以下のように"status"コマンドを実行します．
(@v1.8) pkg> status
Status `~/.julia/environments/v1.8/Project.toml`
[336ed68f] CSV v0.10.4
[861a8166] Combinatorics v1.0.2
[a93c6f00] DataFrames v1.3.5
⌅ [31c24e10] Distributions v0.23.4
[c91e804a] Gadfly v1.3.4
[7073ff75] IJulia v1.23.3
[4076af6c] JuMP v1.3.0
[5424a776] Mamba v0.12.5
[429524aa] Optim v1.7.2
[91a5bcdd] Plots v1.32.1
[438e738f] PyCall v1.94.1
[d330b81b] PyPlot v2.11.0
[ce6b1742] RDatasets v0.7.7
[2913bbd2] StatsBase v0.33.21
⌃ [f3b207a7] StatsPlots v0.15.1
[4b9e565b] SumOfSquares v0.6.2
[37e2e46d] LinearAlgebra
[10745b16] Statistics
Info Packages marked with ⌃ and ⌅ have new versions available, but those with ⌅ cannot be upgraded. To see why use `status --outdated`
2022年7月26日火曜日
ベータ分布
Pythonによるベイズ統計学入門に掲載されているコードを実行してみた際のメモです．
ベータ分布を描画するPythonのコード．
# -*- coding: utf-8 -*-
#%% NumPyの読み込み
import numpy as np
# SciPyのstatsモジュールの読み込み
import scipy.stats as st
# MatplotlibのPyplotモジュールの読み込み
import matplotlib.pyplot as plt
#%% ベータ分布の確率密度関数
q = np.linspace(0, 1, 250)
value_a = np.array([0.5, 1.0, 2.0, 4.0])
value_b = np.array([0.5, 1.0, 2.0, 4.0])
rows = value_a.shape[0]
cols = value_b.shape[0]
fig, ax = plt.subplots(rows, cols, sharex='all', sharey='all',
num=1, facecolor='w')
ax[0, 0].set_xlim(0.0, 1.0)
ax[0, 0].set_ylim(0.0, 4.5)
for row_index in range(rows):
a = value_a[row_index]
ax[row_index, 0].set_ylabel('$\\alpha$ = {0:3.1f}'.format(a),
fontsize=12)
for column_index in range(cols):
b = value_b[column_index]
ax[row_index, column_index].plot(q, st.beta.pdf(q, a, b), 'k-')
if row_index == 0:
ax[0, column_index].set_title('$\\beta$ = {0:3.1f}'.format(b),
fontsize=12)
plt.tight_layout()
plt.savefig('pybayes_fig_beta_distribution.png', dpi=300)
plt.show()
このコードを実行すると，以下の図が描画されます．この図はα，βを変えた時のベータ分布の形状です．
2022年7月25日月曜日
ベルヌーイ分布の成功確率の事前分布（一様分布，ベルヌーイ分布）
Pythonによるベイズ統計学入門に掲載されているコードを実行してみた際のメモです．
ベルヌーイ分布の成功確率の事前分布（一様分布とベルヌーイ分布）を描画するPythonのコード．
# -*- coding: utf-8 -*-
#%% NumPyの読み込み
import numpy as np
# SciPyのstatsモジュールの読み込み
import scipy.stats as st
# MatplotlibのPyplotモジュールの読み込み
import matplotlib.pyplot as plt
# 日本語フォントの設定
from matplotlib.font_manager import FontProperties
import sys
if sys.platform.startswith('win'):
FontPath = 'C:\\Windows\\Fonts\\meiryo.ttc'
elif sys.platform.startswith('darwin'):
FontPath = '/System/Library/Fonts/ヒラギノ角ゴシック W4.ttc'
elif sys.platform.startswith('linux'):
FontPath = '/usr/share/fonts/truetype/takao-gothic/TakaoPGothic.ttf'
else:
print('このPythonコードが対応していないOSを使用しています．')
sys.exit()
jpfont = FontProperties(fname=FontPath)
#%% ベルヌーイ分布の成功確率qの事前分布
fig1 = plt.figure(num=1, facecolor='w')
q = np.linspace(0, 1, 250)
plt.plot(q, st.uniform.pdf(q), 'k-')
plt.plot(q, st.beta.pdf(q, 4, 6), 'k--')
plt.xlim(0, 1)
plt.ylim(0, 2.8)
plt.legend(['(A) 一様分布 ($\\alpha$ = 1, $\\beta$ = 1)',
'(B) ベータ分布 ($\\alpha$ = 4, $\\beta$ = 6)'],
loc='best', frameon=False, prop=jpfont)
plt.xlabel('成功確率 q', fontproperties=jpfont)
plt.ylabel('確率密度', fontproperties=jpfont)
plt.savefig('pybayes_fig_beta_prior.png', dpi=300)
plt.show()
これを実行すると，以下の図が表示されます．
コードの簡単な解説は以下の通りです．
# -*- coding: utf-8 -*-
全角文字である日本語を使用する場合には，様々な文字コードがある（JIS，EUC-JP，Shift JISなど）ので，どの文字コードを使用して作成したコードなのかを明記しておきます．
import numpy as np
# SciPyのstatsモジュールの読み込み
import scipy.stats as st
# MatplotlibのPyplotモジュールの読み込み
import matplotlib.pyplot as pot
NumPy, SciPy, MatPlotLibを読み込んでいます．
# 日本語フォントの設定
from matplotlib.font_manager import FontProperties
import sys
if sys.platform.startswith('win'):
FontPath = 'C:\\Windows\\Fonts\\meiryo.ttc'
elif sys.platform.startswith('darwin'):
FontPath = '/System/Library/Fonts/ヒラギノ角ゴシック W4.ttc'
elif sys.platform.startswith('linux'):
FontPath = '/usr/share/fonts/truetype/takao-gothic/TakaoPGothic.ttf'
else:
print('このPythonコードが対応していないOSを使用しています．')
sys.exit()
jpfont = FontProperties(fname=FontPath)
表示する日本語フォントを選択するために，
import sys
でパッケージ sys を読み込んで，使用しているシステムのOSを判定するための関数
sys.platform.startswith()
を使えるようにしています．この関数は，システムのOSの名前がある文字列（e.g. Windowsならばwin）で始まっていれば真（True），そうでなければ偽（False）を返します．
これを利用して，システムのOSを判定して，OSに標準で付属する日本語フォントをpyplotの作図関数が認識できるようにする命令です．
#%% ベルヌーイ分布の成功確率qの事前分布
fig1 = plt.figure(num=1, facecolor='w')
作図を行うための空白のキャンパスを用意しています．
plt.figure(num=1, facecolor='w)
のnum=1は図の番号を1にするという意味です．特に指定しない場合はpyplot関数です．
plt.figure()
を呼び出すために自動的に通し番号が振られます．
次の facecolor = 'w' というオプションは白地（w: white）のグラフを作成するためのものです．
q = np.linspace(0, 1, 250)
一様分布とベータ分布のグラフを描くための成功確率 q の値を0から1まで等間隔に変化させたグリッドを生成し，それを1次元のNumPy配列qに格納しています．
plt.plot(q, st.uniform.pdf(q), 'k-')
plt.plot(q, st.beta.pdf(q, 4, 6), 'k--')
関数 plt.plot()は，いくつかの点を2次元座標上に打ち，各点を線で結んでグラフにします．なので，各点の横軸と縦軸の座標データとしてplt.plot()に与えなければなりません．これは1次元のNumPy配列の形でplt.plot()に渡されます．
上の行では，一様分布のグラフ（uniform）を作成しています．
先に読み込んでいる
import scipy.stats as at
には，一様分布の確率密度関数 st.uniform.pdf() が含まれているので，これを利用して一様分布の確率密度を計算します．描画の際には'k-'で黒の線を指定しています．
次の行では，ベータ分布の確率密度を計算するため stats 関数 st.beta.pdf() を使用しています．描画の際には'k--'で黒の波線を指定しています．
plt.xlim(0, 1)
plt.ylim(0, 2.8)
関数plt.xlim()は，横軸（x軸）の範囲を指定する関数で，最初の数値が下限，次の数値が上限です．横軸はベルヌーイ分布の成功確率qなので，0を下限，1を上限としています．
plt.ylim() は縦軸（y軸）の範囲を指定する関数で，最初の数値が下限，次の数値が上限です．
plt.legend(['(A) 一様分布 ($\\alpha$ = 1, $\\beta$ = 1)',
'(B) ベータ分布 ($\\alpha$ = 4, $\\beta$ = 6)'],
loc='best', frameon=False, prop=jpfont)
関数plt.legend()は与えられた文字列を要素とするPythonのリスト（これはNumPy配列ではありません）を使って凡例を作成する関数です．
凡例をキレイに表示するために，3つのオプションがあり，'best'は最も見やすい場所に判例を自動的に配置する指示です．明示的に配置する場所を決める際には，'upper left'，upper right'，lower left' などとlocオプションを変えることで変更可能です．
frameon=False は凡例を囲む枠を省くというオプションです．デフォルトでは，凡例を囲む枠が表示されるので，Falseとして枠をなくしています．
prop=jpfontは，日本語フォントを凡例の表示に使用するためのオプションです．
plt.xlabel('成功確率 q', fontproperties=jpfont)
plt.ylabel('確率密度', fontproperties=jpfont)
関数 plt.xlabel() は横軸のラベルを指定する関数で，関数 plt.ylabel()は縦軸のラベルを指定する関数です．
ここでも，ラベルに使われる日本語フォントとして fontproperties = jpfont というオプションで指定しています．
plt.savefig('pybayes_fig_beta_prior.png', dpi=300)
plt.show()
作成した図を画像ファイルに保存したい際には，plt.savefig()の 'pybayes.fig.beta_prior.png' は保存する図の名前，形式（.png）です．
登録:
投稿
(
Atom
)