ガウス分布と指数分布の混合分布

1次元のガウス分布と指数分布の混合分布は以下の式で表される

ここに潜在変数$z_i$を導入すると

と積の形でかけ、尤度計算がしやすい

混合分布の事後分布

事前分布を以下のように定める

このとき、変分ベイズの事後分布は

また、各パラメータについて、以下の式が成り立つ。

各期待値は以下のように表せる。

変分ベイズ学習を実践

$λ=1$の指数分布と$μ=5$のガウス分布を1:1の割合で混合した分布をpythonで発生させる（ガウス分布の分散既知とする）。

この分布を変分ベイズ推定すると以下のような結果が得られる。

左から順に
$λ$、　$μ$、指数分布からサンプリングされる割合
となっており分布のパラメータを学習できていることが窺える

プログラム

import numpy as np
import math
import matplotlib.pyplot as plt
import scipy.special as sp
def F_Ni(zi_rz):
    s = 0
    for z in zi_rz:
        s += z
    return s
def F_xi(Ni,zi_rz,x):
    s = 0
    for i in range(len(x)):
        s += zi_rz[i]*x[i]
    return s
def F_z1(x,E_lambda,E_log_alpha_1,E_log_lambda):
    N = len(x)
    res = [0 for i in range(N)]
    for i in range(N):
        res[i] = math.exp(-x[i]*E_lambda+E_log_alpha_1+E_log_lambda)
    return res
def F_z2(gamma_const,E_mom2,E_log_alpha_2):
    N = len(E_mom2)
    res = [0 for i in range(N)]
    for i in range(N):
        res[i] = math.exp(-gamma_const*E_mom2[i]+E_log_alpha_2)
    return res
def F_z1_rz(z1,z2):
    N = len(z1)
    res = [0 for i in range(N)]
    for i in range(N):
        res[i] = z1[i]/(z1[i]+z2[i])
    return res
def F_z2_rz(z1,z2):
    N = len(z1)
    res = [0 for i in range(N)]
    for i in range(N):
        res[i] = z2[i]/(z1[i]+z2[i])
    return res
def F_E_log_alpha_1(phai1_after_tmp,phai2_after_tmp,N1,N2):
    return sp.digamma(phai1_after_tmp+N1)-sp.digamma(phai1_after_tmp+phai2_after_tmp+N1+N2)
def F_E_log_alpha_2(phai1_after_tmp,phai2_after_tmp,N1,N2):
    return sp.digamma(phai2_after_tmp+N2)-sp.digamma(phai1_after_tmp+phai2_after_tmp+N1+N2)
def F_E_mom2(gamma_after,x,mu_after):
    N = len(x)
    res = [0 for i in range(N)]
    for i in range(N):
        res[i] = 1/(2*gamma_after) + (x[i]-mu_after)**2
    return res
def F_E_lambda(k_after,l_after):
    return k_after/l_after
def F_E_log_lambda(k_after,l_after):
    return sp.digamma(k_after)-math.log(l_after)
#ベイズ更新
def F_k_after(N1,k):
    return N1+k
def F_l_after(N1,x1,l):
    return x1+l
def F_mu0_after(gamma0,mu0,gamma_const,N2,x2):
    return (gamma0*mu0+gamma_const*x2)/(gamma0+gamma_const*N2)
def F_gamma0_after(gamma0,gamma_const,N2):
    return gamma0+gamma_const*N2
def solver(batch,x,params,gamma_const):
    #params = [k_after,l_after,mu0_after,gamma0_after]
    k_after = params[0]
    l_after = params[1]
    mu0_after = params[2]
    gamma0_after = params[3]
    phai1_after = params[4]
    phai2_after = params[5]
    k_after_tmp = params[0]
    l_after_tmp = params[1]
    mu0_after_tmp = params[2]
    gamma0_after_tmp = params[3]
    phai1_after_tmp = params[4]
    phai2_after_tmp = params[5]
    E_lambda = 1
    E_log_alpha_1 = -1
    E_log_alpha_2 = -1
    E_log_lambda = 0
    E_lambda_before = 0
    E_log_alpha_1_before = -1
    E_log_alpha_2_before = -1
    E_log_lambda_before = 0
    z1 = [ 0 for i in range(batch)]
    z2 = [ 0 for i in range(batch)]
    E_mom2 = [ 1 for i in range(batch)]
    E_mom2_before = [ 0 for i in range(batch)]
    eps = 0.00001
    f1 = True
    f2 = True
    f3 = True
    f4 = True
    f5 = True
    while f1 or f2 or f3 or f4 or f5:
        f5 = False
        z1 = F_z1(x,E_lambda,E_log_alpha_1,E_log_lambda)
        z2 = F_z2(gamma_const,E_mom2,E_log_alpha_2)
        z1_rz = F_z1_rz(z1,z2)
        z2_rz = F_z2_rz(z1,z2)
        N1 = F_Ni(z1_rz)
        N2 = F_Ni(z2_rz)
        x1 = F_xi(N1,z1_rz,x)
        x2 = F_xi(N2,z2_rz,x)
        E_log_alpha_1 = F_E_log_alpha_1(phai1_after_tmp,phai2_after_tmp,N1,N2)
        E_log_alpha_2 = F_E_log_alpha_2(phai1_after_tmp,phai2_after_tmp,N1,N2)
        k_after_tmp = F_k_after(N1,k_after_tmp)
        l_after_tmp = F_l_after(N1,x1,l_after_tmp)
        mu0_after_tmp = F_mu0_after(gamma0_after_tmp,mu0_after_tmp,gamma_const,N2,x2)
        gamma0_after_tmp = F_gamma0_after(gamma0_after_tmp,gamma_const,N2)
        E_mom2 = F_E_mom2(gamma0_after_tmp,x,mu0_after_tmp)
        #print(N2)
        #input()
        E_lambda = F_E_lambda(k_after_tmp,l_after_tmp)
        E_log_lambda = F_E_log_lambda(k_after_tmp,l_after_tmp)
        f1 = (abs(E_lambda - E_lambda_before) > eps)
        f2 = (abs(E_log_alpha_1 - E_log_alpha_1_before) > eps)
        f3 = (abs(E_log_alpha_2 - E_log_alpha_2_before) > eps)
        f4 = (abs(E_log_lambda - E_log_lambda_before) > eps)
        #print(f1,f4,f5)
        for i in range(batch):
            f5 += (abs(E_mom2[i] - E_mom2_before[i]) > eps)
        E_lambda_before = E_lambda
        E_log_alpha_1_before = E_log_alpha_1
        E_log_alpha_2_before = E_log_alpha_2
        E_log_lambda_before = E_log_lambda
        for i in range(batch):
            E_mom2_before[i] = E_mom2[i]
    k_after = F_k_after(N1,k_after)
    l_after = F_l_after(N1,x1,l_after)
    mu0_after = F_mu0_after(gamma0_after,mu0_after,gamma_const,N2,x2)
    gamma0_after = F_gamma0_after(gamma0_after,gamma_const,N2)
    return [k_after,l_after,mu0_after,gamma0_after,N1+phai1_after,N2+phai2_after]
def execute_analysis(mixed_samples):
    k_after = 1
    l_after = 1
    mu0_after = 0
    gamma0_after = 1
    gamma_const = 2
    phai1 = 1
    phai2 = 1
    params = [k_after,l_after,mu0_after,gamma0_after,phai1,phai2]
    batch = 10
    x = [ 0 for i in range(batch)]
    for i in range(int(len(mixed_samples)/batch)):
        x = [mixed_samples[batch*i+j] for j in range(batch)]
        params = solver(batch,x,params,gamma_const)
        #params[0] /= 100
        #params[1] /= 100
        print(params[0]/params[1],params[2],(params[4]-1)/(params[4]+params[5]-2))
# パラメータ設定
mu = 5  # ガウス分布の平均
sigma = 0.5  # ガウス分布の標準偏差
lambda_exp = 1  # 指数分布のパラメータ (レート λ)
# ガウス分布と指数分布の混合割合 (0.5はそれぞれ半分ずつ)
mixing_ratio = 0.5
# 乱数の個数
num_samples = 10000
# どちらの分布からサンプルするかを決めるフラグを生成 (0: ガウス, 1: 指数)
flags = np.random.choice([0, 1], size=num_samples, p=[mixing_ratio, 1 - mixing_ratio])
# ガウス分布からのサンプル
gaussian_samples = np.random.normal(mu, sigma, num_samples)
# 指数分布からのサンプル
exponential_samples = np.random.exponential(1 / lambda_exp, num_samples)
# フラグに基づいて最終的なサンプルを決定
mixed_samples = np.where(flags == 0, gaussian_samples, exponential_samples)
# ヒストグラムを描画
plt.figure(figsize=(8, 6))
plt.hist(mixed_samples, bins=30, alpha=0.7, color='b', edgecolor='black')
plt.title('Mixed Distribution: Gaussian and Exponential')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.grid(True)
# グラフ表示
plt.show()
execute_analysis(mixed_samples)