<jupyter_start><jupyter_text>Hyperbolic Secant Mixture modelの推定+ モデル:    + $p(x|w) = \sum_{k=1}^K a_k \frac{\sqrt{s_k}}{2} \frac{1}{\cosh(\frac{s_k}{2}(x - b_k))}$    + $x, b_k \in \mathbb{R}, s_k \in \mathbb{R}_+$+ 事前分布:    + $\varphi(w) = Dir(a|\{ \alpha_k \}_{k=1}^K) \prod_{k=1}^K N(b_k|0, (s_k \beta_k)^{-1} ) Gam(s_k|\gamma_k, \delta_k)$+ アルゴリズム   1. 以下の初期化を行う:       + $g_{ik}(\eta), v_{ik}(\eta), h_{ik}(\xi), u_{ik}(\xi)$   2. 以下を計算する       + $\hat{\alpha}_k = \alpha_k + \sum_{i=1}^n u_{ik}(\xi)$       + $\hat{\beta}_k = \beta_k + \sum_{i=1}^n -2v_{ik}(\eta)$       + $\hat{m}_k = \frac{1}{\hat{\beta}_k} \sum_{i=1}^n -2v_{ik}(\eta)x_i$       + $\hat{\gamma}_k = \gamma_k + \frac{1}{2}\sum_{i=1}^n u_{ik}(\xi)$       + $\hat{\delta}_k = \delta_k + \sum_{i=1}^n -v_{ik}(\eta)x_i^2 - \frac{\hat{\beta}_k}{2}\hat{m}_k^2$   3. 以下を計算する       + $g_{ik}(\eta) = \frac{\hat{\gamma}_k}{\hat{\delta}_k} (x_i - \hat{m}_k)^2 + \frac{1}{\hat{\beta}_k}$       + $v_{ik}(\eta) = -u_{ik}(\xi)\frac{ \tanh(\sqrt{g_{ik}(\eta)}/2) }{4\sqrt{g_{ik}(\eta)}}$   4. 以下を計算する       + $h_{ik}(\xi) = \psi(\hat{\alpha}_k) - \psi(\sum_{l=1}^K \hat{\alpha}_l) + \frac{1}{2}(\psi(\hat{\gamma}_k) - \log(\hat{\delta}_k)) - \log(\cosh(\sqrt{g_{ik}(\eta)}/2)) $       + $u_{ik}(\xi) = \frac{ \exp(h_{ik}(\xi)) }{ \sum_{l=1}^K \exp(h_{il}(\xi)) }$       + ここで,$\psi(x) = \frac{d}{dx}\log \Gamma(x)$   5. 2.に戻る   + 評価関数$\overline{F}_{\xi, \eta}(x^n)$の値    + $\overline{F}_{\xi, \eta}(x^n) = - \phi(h(\xi)) - \psi(g(\eta)) + u(\xi) \cdot h(\xi) + v(\eta) \cdot g(\eta) $          $+ n \log 2 \pi + \log \Gamma(\sum_{l = 1}^K \hat{\alpha}_l) - \log \Gamma({\sum_{l = 1}^K\alpha}_l) + \sum_{k=1}^K \log \frac{\Gamma(\alpha_k)}{\Gamma(\hat{\alpha}_k)}$          $+ \sum_{k=1}^K \bigl\{ \frac{1}{2} \log \frac{\hat{\beta}_k}{\beta_k} + \hat{\gamma}_k \log \hat{\delta}_k - \gamma_k \log \delta_k - \log \Gamma(\hat{\gamma}_k) + \log \Gamma(\gamma_k) \bigr\}$<jupyter_code>%matplotlib inline
import numpy as np
from scipy.special import gammaln, psi
import matplotlib.pyplot as plt
import seaborn as sns<jupyter_output><empty_output><jupyter_text>HSMMに従う確率変数を生成する<jupyter_code>data_seed = 20190505
n = 1000
true_ratio = np.array([0.5, 0.5])
true_s = np.array([4, 0.5])
true_b = np.array([4, -4])
np.random.seed(data_seed)
true_label = np.random.multinomial(n = 1, pvals = true_ratio, size = n)
true_label_arg = np.argmax(true_label, axis = 1)
def random_hsm(n, loc = 0, scale = 1):
    """
    hyperbolic secant distributionに従う確率変数を生成する
    $Y \sim standard_cauchy(x)$の時,  
    $X = \frac{2}{s}\sinh^{-1}(Y) + b$は,  
    $X \sim p(x) = \frac{s}{2\pi}\frac{1}{\cosh(s(x-b)/2)}$に従うはずなので、その確認を行う
    """
    Y = np.random.standard_cauchy(size=n)
    X = 2/np.sqrt(scale)*np.arcsinh(Y) + loc    
    return X
x = np.zeros(n)
for i in range(n):
    x[i] = random_hsm(n=1, loc=true_b[true_label_arg[i]], scale=true_s[true_label_arg[i]])
plt.scatter(x,np.zeros(n))
plt.hist(x)
plt.hist(x[np.where(true_label_arg == 1)])
plt.hist(x[np.where(true_label_arg == 0)])<jupyter_output><empty_output><jupyter_text>学習を行う<jupyter_code>### 初期設定を行う
iteration = 1000
learning_seed = 20190510

## 学習モデルのコンポーネントの数
K = 3

## ハイパーパラメータ
pri_alpha = 3
pri_beta = 0.001
pri_gamma = 2
pri_delta = 2

np.random.seed(learning_seed)

energy = np.zeros(iteration)
### 初期値の設定
est_h_xi = np.random.gamma(shape = 1, size = (n, K))
est_u_xi = np.exp(est_h_xi) / np.repeat(np.exp(est_h_xi).sum(axis = 1), K).reshape(n,K)
# est_u_xi = np.random.dirichlet(alpha = np.ones(K), size=n)

est_g_eta = np.abs(np.random.normal(size=(n,K)))
est_v_eta = - est_u_xi * np.tanh(np.sqrt(est_g_eta)/2)/(4*np.sqrt(est_g_eta))
### 繰り返し計算部分
for ite in range(iteration):
    ### 事後分布の更新
    est_alpha = pri_alpha + est_u_xi.sum(axis = 0)
    est_beta = pri_beta + (-2*est_v_eta.sum(axis = 0))
    est_m = (-2*x @ est_v_eta)/est_beta
    est_gamma = pri_gamma + est_u_xi.sum(axis = 0)/2
    est_delta = pri_delta - (x**2) @ est_v_eta - est_beta / 2 * est_m**2
    ### 対数双曲線関数の更新
    est_g_eta = est_gamma / est_delta * (np.repeat(x,K).reshape(n,K) - np.repeat(est_m,n).reshape(K,n).T)**2 + 1/est_beta
    est_v_eta = - est_u_xi * np.tanh(np.sqrt(est_g_eta)/2)/(4*np.sqrt(est_g_eta))
    
    ### 隠れ変数の更新
    est_h_xi = np.repeat(psi(est_alpha) - psi(est_alpha.sum()) + (psi(est_gamma) - np.log(est_delta))/2, n).reshape(K,n).T - np.log(np.cosh(np.sqrt(est_g_eta)/2))
    max_h_xi = est_h_xi.max(axis = 1)
    norm_h_xi = est_h_xi - np.repeat(max_h_xi,K).reshape(n,K)
    est_u_xi = np.exp(norm_h_xi) / np.repeat(np.exp(norm_h_xi).sum(axis = 1), K).reshape(n,K)

    ### エネルギーの計算
    energy[ite] = (est_u_xi * np.log(np.cosh(np.sqrt(est_g_eta)/2))).sum() - (np.log(np.exp(norm_h_xi).sum(axis = 1)) + max_h_xi).sum() + (est_u_xi * est_h_xi + est_v_eta * est_g_eta).sum()
    energy[ite] += gammaln(est_alpha.sum()) - gammaln(K*pri_alpha) + (-gammaln(est_alpha) + gammaln(pri_alpha)).sum()
    energy[ite] += (np.log(est_beta/pri_beta)/2 + est_gamma * np.log(est_delta) - pri_gamma * np.log(pri_delta) - gammaln(est_gamma) + gammaln(pri_gamma)).sum()
    print(energy[ite])
{
    "est_ratio": est_alpha / est_alpha.sum(),
    "est_b": est_m,
    "est_s": est_gamma / est_delta
}<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.23710944987540733, 'fraction_numerical': 0.023193406172129578, 'mean_word_length': 3.7522768670309654, 'pattern_counts': {'":': 3, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '102000', 'license': 'MIT', 'n_tokens_mistral': 2635, 'n_tokens_neox': 2562, 'n_words': 347, 'path': 'HSMM/HSMM_test.ipynb', 'repo_name': 'fumish/lva', 'size': '5217'}