Python で ガンマ分布の計算をする

import numpy as np
shape, scale = 2., 2.
s = np.random.gamma(shape, scale, 1000)
print(s)

%matplotlib inline
import numpy as np
shape, scale = 1., 1
s = np.random.gamma(shape, scale, 10000)
import matplotlib.pyplot as plt
import scipy.special as sps
count, bins, ignored = plt.hist(s, 50, density=True)
y = bins**(shape-1)*(np.exp(-bins/scale) /
                      (sps.gamma(shape)*scale**shape))
plt.plot(bins, y, linewidth=2, color='r')

%matplotlib inline
import numpy as np
import pylab as plt

target = 1
beta = 1.0/target
Y = np.random.exponential(beta, 10000)
plt.hist(Y, density=True, bins=200,lw=0,alpha=.8)
plt.show()

%matplotlib inline
import matplotlib.pyplot as plt

def plot_gamma(shape, scale):
    s = np.random.gamma(shape, scale, 1000)
    import scipy.special as sps
    count, bins, ignored = plt.hist(s, 50, density=True)
    y = bins**(shape-1)*(np.exp(-bins/scale) /
                      (sps.gamma(shape)*scale**shape))
    plt.plot(bins, y, linewidth=2, color='r')

shape, scale = 1., 1.
plot_gamma(shape, scale)
shape, scale = 2., 1.
plot_gamma(shape, scale)
shape, scale = 3., 1.
plot_gamma(shape, scale)
import matplotlib.pyplot as plt
plt.show()

%matplotlib inline
from scipy.stats import gamma
import matplotlib.pyplot as plt
import numpy as np


def plot_gamma(a):
    fig, ax = plt.subplots(1, 1)
    # -----------------------------
    # モーメントを求める
    # ---------------------
    # 以下をいつも忘れる
    # mean 平均, var 分散,skew 歪度,kurt 尖度

    mean, var, skew, kurt = gamma.stats(a, moments='mvsk')
    print("平均: {0}, 分散: {1}, 歪度: {2}, 尖度: {3}".format(mean, var, skew, kurt))

    # numpy.linespace で等差数列を生成する
    # gamma.ppf は、％点関数
    x = np.linspace(gamma.ppf(0.01, a),
                     gamma.ppf(0.99, a), 100)

    #  pdf は、確率密度関数 これもいつも忘れる。   
    # plot の引数 alpha は、透過度を示す値、lw は線の太さを指定     
    ax.plot(x, gamma.pdf(x, a),'r-', lw=5, alpha=0.6, label='gamma pdf')

    # 固定型統計分布: frozen RV
    rv = gamma(a)
    ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')

    #  xが0.001,0.5, 0.999 での％点関数 を求める    
    vals = gamma.ppf([0.001, 0.5, 0.999], a)

    # numpy.allclose で配列の等価、不等価の判定ができる。    
    allclose_result = np.allclose([0.001, 0.5, 0.999], gamma.cdf(vals, a))
    print("allclose_result: {0}".format(allclose_result))

    # rvs でガンマ分布に従う乱数を生成、size は 1000 なので 要素数 1000 になる。 デフォルトは size = 1
    r = gamma.rvs(a, size=1000)
    # ヒストグラム描画、bins 指定なしでデフォルト 10 , stepfilled は 塗りつぶしありの線    
    # densityがTrueの場合、戻確率密度を形成するように正規化される。
    ax.hist(r, density=True, histtype='stepfilled', alpha=0.2)
    ax.legend(loc='best', frameon=False)
    plt.show()


# a=5 のガンマ分布
plot_gamma(5.0)

# a=1  指数分布
plot_gamma(1.0)

%matplotlib inline
import tensorflow as tf
import pylab as plt
#  10000要素 alpha=1 つまり指数分布に従う配列を作成する     
samples = tf.random_gamma(shape=[10000], alpha=1)

# Sessionオブジェクトを作る
# 初めてTensorFLow 触りましたが、考え方に慣れておらずあまり意味わかっておりません。    
with tf.Session() as sess:
    print(samples.__dict__)

    y = samples.eval()
    plt.hist(y, density=True, bins=200,lw=0,alpha=.8)
    plt.show()

# alpha を定数化
alpha = tf.constant([[1.],[3.]])
# beta を定数化
beta = tf.constant([[3., 4.]])
# 2x2 の多次元配列を、3要素生成する
samples = tf.random_gamma([3], alpha=alpha, beta=beta)
with tf.Session() as sess:    
    y = samples.eval()
    print(y)


# tf.reduce_mean 与えたリストに入っている数値の平均値を求める関数    
# tf.square 要素ごとに二乗をとる
loss = tf.reduce_mean(tf.square(samples))
# 勾配を求める    
dloss_dalpha, dloss_dbeta = tf.gradients(loss, [alpha, beta])
# 損失関数の不偏確率導関数?    
print("alpha.shape == dloss_dalpha.shape", alpha.shape == dloss_dalpha.shape)  # True
print("beta.shape == dloss_dbeta.shape", beta.shape == dloss_dbeta.shape)  # True