ガンマ分布

いつの間にか社会人です。
最近業務でいくつかの確率分布を触っているのですが、パラメータが変わると分布がどう変わるのか、パッと知りたかったので

ガンマ分布を例に、PythonのBokehライブラリを使ってインタラクティブなグラフを作ってみました。
実際に動かせるので触ってみてください。
※PCサイトでないと表示されないみたいです。

ちなみにコードは以下のような感じです。JupyterNotebook上で実行してください。
最後のbokeh.plotting.show()の代わりにbokeh.embed.components()を使ってあげればscriptタグとdivタグを吐き出してくれます。
これに加えてブログのヘッダーにBokeh.jsを読み込むコードを追加すればOKです。詳しくは公式のドキュメントを見てください。

import numpy as np
import scipy.special as special
import bokeh
import bokeh.plotting as plt
import bokeh.models as models
import bokeh.layouts as layouts
import bokeh.embed as embed

print(bokeh.__version__)

def plot_gamma_distribution(x, k=2, beta=1):
    gamma_y = beta**k / special.gamma(k) * x**(k-1) * np.exp(-beta * x)
    power_y = x**(k-1)
    exp_y = np.exp(-beta * x)
    source = plt.ColumnDataSource(data=dict(x=x, gamma=gamma_y, power=power_y, exp=exp_y))
    
    gamma = plt.figure(plot_height=200, plot_width=300, title='ガンマ分布', y_range=(0, 1))
    power = plt.figure(plot_height=200, plot_width=300, title='べき函数部分')
    exp = plt.figure(plot_height=200, plot_width=300, title='指数関数部分')
    gamma.line('x', 'gamma', source=source)
    power.line('x', 'power', source=source)
    exp.line('x', 'exp', source=source)
    
    shape_slider = models.Slider(start=0.1, end=10, value=k, step=0.01, title='shape parameter')
    scale_slider = models.Slider(start=0.1, end=10, value=beta, step=0.01, title='scale parameter')
    
    callback = models.CustomJS(args=dict(source=source, shape=shape_slider, scale=scale_slider), code="""
    function Gamma(x) {
    var p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
        771.32342877765313, -176.61502916214059, 12.507343278686905,
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
    ];
 
    var g = 7;
    if (x < 0.5) {
        return Math.PI / (Math.sin(Math.PI * x) * Gamma(1 - x));
    }
 
    x -= 1;
    var a = p[0];
    var t = x + g + 0.5;
    for (var i = 1; i < p.length; i++) {
        a += p[i] / (x + i);
    }
 
    return Math.sqrt(2 * Math.PI) * Math.pow(t, x + 0.5) * Math.exp(-t) * a;
    }
    
    const data = source.data;
    const k = shape.value;
    const beta = scale.value;
    const x = data['x'];
    const gamma = data['gamma'];
    const power = data['power'];
    const exp = data['exp'];
    for (var i = 0; i < x.length; i++) {
        gamma[i] = beta**k / Gamma(k) * x[i]**(k-1) * Math.exp(-beta * x[i]);
        power[i] = x[i]**(k-1);
        exp[i] = Math.exp(-beta * x[i]);
    }
    source.change.emit();
    """)
    
    shape_slider.js_on_change('value', callback)
    scale_slider.js_on_change('value', callback)
    
    layout = layouts.row(layouts.column(gamma, power, exp), layouts.column(shape_slider, scale_slider))
    return layout

plt.output_notebook()
fig = plot_gamma_distribution(np.linspace(0, 10, 1001))
plt.show(fig)

JSでのガンマ関数の実装は

Gamma function - Rosetta Code

こちらのコードをそのままコピペしました。