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

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

[ 3.82769986  0.57987294  8.71851465  2.09624801  3.07992525  1.33682621
  1.52681352  1.1038436   0.32452196  2.61687864  2.31279027  4.1192877
 11.93373265  1.49622973  2.43055922  1.36888852  2.43573475  6.16885654
  3.43197913  4.41584039  2.63252622  5.84375382  2.4926901   4.09717026
  0.54727079  1.37837228  1.57989827  3.60335118  7.90316212  1.73799822
  3.02557085  1.22738672  0.64397655  6.65220901  7.68344347  9.00825348
  1.30000658  4.21622653  3.81001449  5.44076707  0.29190632  3.4004747
  3.31920514  2.61953315  4.02076256  1.57446158  4.48032433  3.14871272
  0.77085133  8.99744653  1.55685237  1.53318994  4.99302374  1.25386129
  4.28039121  5.74574233  1.36186417 10.10874118  3.65846922  2.02873105
  3.34896627  0.52464966  2.091951    2.28638439  2.75070701 15.00023712
  1.66613084]

%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.lines.Line2D at 0x130fff2b0>]

%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)

平均: 5.0, 分散: 5.0, 歪度: 0.8944271909999159, 尖度: 1.2
allclose_result: True

平均: 1.0, 分散: 1.0, 歪度: 2.0, 尖度: 6.0
allclose_result: True

%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

{'_op': <tf.Operation 'random_gamma_27/Maximum' type=Maximum>, '_value_index': 0, '_dtype': tf.float32, '_shape_val': None, '_consumers': [], '_id': 318}

[[[0.0756161  0.03140004]
  [1.9382031  0.4650259 ]]

 [[0.1251423  0.14494717]
  [1.268082   0.6238722 ]]

 [[0.27819878 0.08588854]
  [0.98764956 0.29059193]]]
alpha.shape == dloss_dalpha.shape True
beta.shape == dloss_dbeta.shape True