updating rule of VB

Author: yuta-nakahara

bayesml/gaussianmixture/_gaussianmixture.py

@@ -8,7 +8,7 @@
 from scipy.stats import wishart as ss_wishart
 from scipy.stats import multivariate_t as ss_multivariate_t
 from scipy.stats import dirichlet as ss_dirichlet
-from scipy.special import gammaln, digamma, xlogy
+from scipy.special import gammaln, digamma, xlogy, logsumexp
 import matplotlib.pyplot as plt
 
 from .. import base
@@ -499,15 +499,28 @@ def __init__(
         self.hn_w_mats = np.empty([self.num_classes,self.degree,self.degree])
         self.hn_w_mats_inv = np.empty([self.num_classes,self.degree,self.degree])
 
-        # statistics
+        self.ln_rho = None
         self.r_vecs = None
-        self.x_bar_vecs = np.empty([self.num_classes,self.degree])
-        self.ns = np.empty(self.num_classes)
-        self.s_mats = np.empty([self.num_classes,self.degree,self.degree])
         self.e_lambda_mats = np.empty([self.num_classes,self.degree,self.degree])
         self.e_ln_lambda_dets = np.empty(self.num_classes)
+        self.ln_b_hn_w_nus = np.empty(self.num_classes)
         self.e_ln_pi_vec = np.empty(self.num_classes)
 
+        # statistics
+        self.x_bar_vecs = np.empty([self.num_classes,self.degree])
+        self.ns = np.empty(self.num_classes)
+        self.s_mats = np.empty([self.num_classes,self.degree,self.degree])
+
+        # variational lower bound
+        self.vl = 0.0
+        self.vl_p_x = 0.0
+        self.vl_p_z = 0.0
+        self.vl_p_pi = 0.0
+        self.vl_p_mu_lambda = 0.0
+        self.vl_q_z = 0.0
+        self.vl_q_pi = 0.0
+        self.vl_q_mu_lambda = 0.0
+
         # p_params
         self.p_pi_vec = np.empty([self.num_classes])
         self.p_mu_vecs = np.empty([self.num_classes,self.degree])
@@ -695,6 +708,21 @@ def reset_hn_params(self):
         self.hn_w_mats[:] = self.h0_w_mats
         self.hn_w_mats_inv = np.linalg.inv(self.hn_w_mats)
 
+        self.e_lambda_mats[:] = self.hn_nus[:,np.newaxis,np.newaxis] * self.hn_w_mats
+        self.e_ln_lambda_dets[:] = (
+            np.sum(digamma((self.hn_nus[:,np.newaxis]-np.arange(self.degree)) / 2.0),axis=1)
+            + self.degree*np.log(2.0)
+            - np.linalg.slogdet(self.hn_w_mats_inv)[1]
+            )
+        self.e_ln_pi_vec[:] = digamma(self.hn_alpha_vec) - digamma(self.hn_alpha_vec.sum())
+        self.ln_b_hn_w_nus[:] = (
+            self.hn_nus*np.linalg.slogdet(self.hn_w_mats_inv)[1]
+            - self.hn_nus*self.degree*np.log(2.0)
+            - self.degree*(self.degree-1)/2.0*np.log(np.pi)
+            - np.sum(gammaln((self.hn_nus[:,np.newaxis]-np.arange(self.degree)) / 2.0),
+                     axis=1) * 2.0
+            ) / 2.0
+
         self.calc_pred_dist()
 
     def overwrite_h0_params(self):
@@ -713,22 +741,8 @@ def overwrite_h0_params(self):
         self.calc_pred_dist()
 
     def calc_vl(self):
-        self.e_lambda_mats = self.hn_nus[:,np.newaxis,np.newaxis] * self.hn_w_mats
-        self.e_ln_lambda_dets = (np.sum(digamma((self.hn_nus[:,np.newaxis]-np.arange(self.degree)) / 2.0),axis=1)
-                            + self.degree*np.log(2.0)
-                            - np.linalg.slogdet(self.hn_w_mats_inv)[1])
-        self.e_ln_pi_vec = digamma(self.hn_alpha_vec) - digamma(self.hn_alpha_vec.sum())
-        
-        # tentative
-        self.ns = np.ones(self.num_classes) * 10
-        self.s_mats = np.tile(np.identity(self.degree),[self.num_classes,1,1]) * 5
-        self.r_vecs = np.ones([20,self.degree])/self.degree
-        self.x_bar_vecs = np.ones([self.num_classes,self.degree])
-
-        vl = 0.0
-
         # E[ln p(X|Z,mu,Lambda)]
-        vl += np.sum(
+        self.vl_p_x = np.sum(
             self.ns
             * (self.e_ln_lambda_dets - self.degree / self.hn_kappas
                - (self.s_mats * self.e_lambda_mats).sum(axis=(1,2))
@@ -741,63 +755,173 @@ def calc_vl(self):
             ) / 2.0
 
         # E[ln p(Z|pi)]
-        vl += (self.ns * self.e_ln_pi_vec).sum()
+        self.vl_p_z = (self.ns * self.e_ln_pi_vec).sum()
 
         # E[ln p(pi)]
-        vl += self.LN_C_H0_ALPHA + ((self.h0_alpha_vec - 1) * self.e_ln_pi_vec).sum()
+        self.vl_p_pi = self.LN_C_H0_ALPHA + ((self.h0_alpha_vec - 1) * self.e_ln_pi_vec).sum()
 
         # E[ln p(mu,Lambda)]
-        vl += np.sum(
-            self.degree * (np.log(self.h0_kappas) - np.log(2*np.pi) - self.h0_kappas/self.hn_kappas)
-            - ((self.hn_m_vecs - self.h0_m_vecs)[:,np.newaxis,:]
-               @ self.e_lambda_mats
-               @ (self.hn_m_vecs - self.h0_m_vecs)[:,:,np.newaxis])[:,0,0]
+        self.vl_p_mu_lambda = np.sum(
+            self.degree * (np.log(self.h0_kappas) - np.log(2*np.pi)
+                           - self.h0_kappas/self.hn_kappas)
+            - self.h0_kappas * ((self.hn_m_vecs - self.h0_m_vecs)[:,np.newaxis,:]
+                                @ self.e_lambda_mats
+                                @ (self.hn_m_vecs - self.h0_m_vecs)[:,:,np.newaxis])[:,0,0]
             + 2.0 * self.LN_B_H0_W_NUS
-            + (self.h0_nus - self.degree) / 2.0 * self.e_ln_lambda_dets
-            - np.sum(self.h0_w_mats_inv * self.hn_w_mats,axis=(1,2))
+            + (self.h0_nus - self.degree) * self.e_ln_lambda_dets
+            - np.sum(self.h0_w_mats_inv * self.e_lambda_mats,axis=(1,2))
             ) / 2.0
 
         # E[ln q(Z|pi)]
-        vl -= np.sum(xlogy(self.r_vecs,self.r_vecs))
+        self.vl_q_z = -np.sum(xlogy(self.r_vecs,self.r_vecs))
 
         # E[ln q(pi)]
-        vl += ss_dirichlet.entropy(self.hn_alpha_vec)
+        self.vl_q_pi = ss_dirichlet.entropy(self.hn_alpha_vec)
 
         # E[ln q(mu,Lambda)]
-        vl +=  np.sum(
+        self.vl_q_mu_lambda =  np.sum(
             + self.degree * (1.0 + np.log(2.0*np.pi) - np.log(self.hn_kappas))
-            - self.LN_B_H0_W_NUS * 2.0
+            - self.ln_b_hn_w_nus * 2.0
             - (self.hn_nus-self.degree)*self.e_ln_lambda_dets
             + self.hn_nus * self.degree
             ) / 2.0
 
-        return vl
+        # print(self.vl_p_x,
+        #       self.vl_p_z,
+        #       self.vl_p_pi,
+        #       self.vl_p_mu_lambda,
+        #       self.vl_q_z,
+        #       self.vl_q_pi,
+        #       self.vl_q_mu_lambda,
+        #       )
+
+        self.vl = (self.vl_p_x
+                   + self.vl_p_z
+                   + self.vl_p_pi
+                   + self.vl_p_mu_lambda
+                   + self.vl_q_z
+                   + self.vl_q_pi
+                   + self.vl_q_mu_lambda)
+
+    def _calc_statistics(self,x):
+        self.ns[:] = self.r_vecs.sum(axis=0)
+        self.x_bar_vecs[:] = (self.r_vecs[:,:,np.newaxis] * x[:,np.newaxis,:]).sum(axis=0) / self.ns[:,np.newaxis]
+        self.s_mats[:] = np.sum(self.r_vecs[:,:,np.newaxis,np.newaxis]
+                                * ((x[:,np.newaxis,:] - self.x_bar_vecs)[:,:,:,np.newaxis]
+                                   @ (x[:,np.newaxis,:] - self.x_bar_vecs)[:,:,np.newaxis,:]),
+                                axis=0) / self.ns[:,np.newaxis,np.newaxis]
+
+    def _init_q_z(self):
+        self.r_vecs[:] = self.rng.dirichlet(np.ones(self.num_classes),self.r_vecs.shape[0])
+
+    def _update_q_pi(self):
+        self.hn_alpha_vec[:] = self.h0_alpha_vec + self.ns
+        self.e_ln_pi_vec[:] = digamma(self.hn_alpha_vec) - digamma(self.hn_alpha_vec.sum())
+
+    def _update_q_mu_lambda(self):
+        self.hn_kappas[:] = self.h0_kappas + self.ns
+        self.hn_m_vecs[:] = (self.h0_kappas[:,np.newaxis] * self.h0_m_vecs
+                             + self.ns[:,np.newaxis] * self.x_bar_vecs) / self.hn_kappas[:,np.newaxis]
+        self.hn_nus[:] = self.h0_nus + self.ns
+        self.hn_w_mats_inv[:] = (self.h0_w_mats_inv
+                                 + self.ns[:,np.newaxis,np.newaxis] * self.s_mats
+                                 + (self.h0_kappas * self.ns / self.hn_kappas)[:,np.newaxis,np.newaxis]
+                                   * ((self.x_bar_vecs - self.h0_m_vecs)[:,:,np.newaxis]
+                                      @ (self.x_bar_vecs - self.h0_m_vecs)[:,np.newaxis,:])
+                                 )
+        self.hn_w_mats[:] = np.linalg.inv(self.hn_w_mats_inv)
+        self.e_lambda_mats[:] = self.hn_nus[:,np.newaxis,np.newaxis] * self.hn_w_mats
+        self.e_ln_lambda_dets[:] = (np.sum(digamma((self.hn_nus[:,np.newaxis]-np.arange(self.degree)) / 2.0),axis=1)
+                            + self.degree*np.log(2.0)
+                            - np.linalg.slogdet(self.hn_w_mats_inv)[1])
+        self.ln_b_hn_w_nus[:] = (
+            self.hn_nus*np.linalg.slogdet(self.hn_w_mats_inv)[1]
+            - self.hn_nus*self.degree*np.log(2.0)
+            - self.degree*(self.degree-1)/2.0*np.log(np.pi)
+            - np.sum(gammaln((self.hn_nus[:,np.newaxis]-np.arange(self.degree)) / 2.0),
+                     axis=1) * 2.0
+            ) / 2.0
 
-    def update_posterior(self,x):
-        pass
-#         """Update the hyperparameters of the posterior distribution using traning data.
+    def _update_q_z(self,x):
+        self.ln_rho[:] = (self.e_ln_pi_vec
+                          + (self.e_ln_lambda_dets
+                             - self.degree * np.log(2*np.pi)
+                             - self.degree / self.hn_kappas
+                             - ((x[:,np.newaxis,:]-self.hn_m_vecs)[:,:,np.newaxis,:]
+                                @ self.e_lambda_mats
+                                @ (x[:,np.newaxis,:]-self.hn_m_vecs)[:,:,:,np.newaxis]
+                                )[:,:,0,0]
+                             ) / 2.0
+                          )
+        self.r_vecs[:] = np.exp(self.ln_rho - logsumexp(self.ln_rho,axis=1,keepdims=True))
+        # self.r_vecs[:] = np.exp(self.ln_rho - self.ln_rho.max(axis=1,keepdims=True))
+        # self.r_vecs[:] /= self.r_vecs.sum(axis=1,keepdims=True)
+
+    def update_posterior(self,x,max_itr=100,num_init=10,tolerance=1.0E-8):
+        """Update the hyperparameters of the posterior distribution using traning data.
 
-#         Parameters
-#         ----------
-#         x : numpy.ndarray
-#             All the elements must be real number.
-#         """
-#         _check.float_vecs(x,'x',DataFormatError)
-#         if self.degree > 1 and x.shape[-1] != self.degree:
-#             raise(DataFormatError(f"x.shape[-1] must be degree:{self.degree}"))
-#         x = x.reshape(-1,self.degree)
-
-#         n = x.shape[0]
-#         x_bar = x.sum(axis=0)/n
-
-#         self.hn_w_mat_inv[:] = (self.hn_w_mat_inv + (x-x_bar).T @ (x-x_bar)
-#                                 + (x_bar - self.hn_m_vec)[:,np.newaxis] @ (x_bar - self.hn_m_vec)[np.newaxis,:]
-#                                   * self.hn_kappa * n / (self.hn_kappa + n))
-#         self.hn_m_vec[:] = (self.hn_kappa*self.hn_m_vec + n*x_bar) / (self.hn_kappa+n)
-#         self.hn_kappa += n
-#         self.hn_nu += n
-
-#         self.hn_w_mat[:] = np.linalg.inv(self.hn_w_mat_inv) 
+        Parameters
+        ----------
+        x : numpy.ndarray
+            All the elements must be real number.
+        max_itr : int, optional
+            maximum number of iterations, by default 100
+        num_init : int, optional
+            number of initializations, by default 10
+        tolerance : float, optional
+            convergence croterion of variational lower bound, by default 1.0E-8
+        """
+        _check.float_vecs(x,'x',DataFormatError)
+        if self.degree > 1 and x.shape[-1] != self.degree:
+            raise(DataFormatError(
+                "x.shape[-1] must be self.degree: "
+                + f"x.shape[-1]={x.shape[-1]}, self.degree={self.degree}"))
+        x = x.reshape(-1,self.degree)
+        self.ln_rho = np.empty([x.shape[0],self.num_classes])
+        self.r_vecs = np.empty([x.shape[0],self.num_classes])
+
+        tmp_vl = 0.0
+        tmp_alpha_vec = np.copy(self.hn_alpha_vec)
+        tmp_m_vecs = np.copy(self.hn_m_vecs)
+        tmp_kappas = np.copy(self.hn_kappas)
+        tmp_nus = np.copy(self.hn_nus)
+        tmp_w_mats = np.copy(self.hn_w_mats)
+        tmp_w_mats_inv = np.copy(self.hn_w_mats_inv)
+
+        for i in range(num_init):
+            self._init_q_z()
+            self._calc_statistics(x)
+            self.calc_vl()
+            print(f'\r{i}. VL: {self.vl}',end='')
+            for t in range(max_itr):
+                vl_before = self.vl
+
+                self._update_q_mu_lambda()
+                self._update_q_pi()
+                self._update_q_z(x)
+                self._calc_statistics(x)
+                self.calc_vl()
+                print(f'\r{i}. VL: {self.vl}',end='')
+                if np.abs((self.vl-vl_before)/vl_before) < tolerance:
+                    break
+            if i==0 or self.vl > tmp_vl:
+                print('*')
+                tmp_vl = self.vl
+                tmp_alpha_vec[:] = self.hn_alpha_vec
+                tmp_m_vecs[:] = self.hn_m_vecs
+                tmp_kappas[:] = self.hn_kappas
+                tmp_nus[:] = self.hn_nus
+                tmp_w_mats[:] = self.hn_w_mats
+                tmp_w_mats_inv[:] = self.hn_w_mats_inv
+            else:
+                print('')
+        
+        self.hn_alpha_vec[:] = tmp_alpha_vec
+        self.hn_m_vecs[:] = tmp_m_vecs
+        self.hn_kappas[:] = tmp_kappas
+        self.hn_nus[:] = tmp_nus
+        self.hn_w_mats[:] = tmp_w_mats
+        self.hn_w_mats_inv[:] = tmp_w_mats_inv        
 
     def estimate_params(self,loss="squared"):
         pass

bayesml/gaussianmixture/gaussianmixture.md

@@ -79,9 +79,9 @@ $$
 \begin{align}
     N_k^{(t)} &= \sum_{i=1}^n r_{i,k}^{(t)}, \\
     \bar{\boldsymbol{x}}_k^{(t)} &= \frac{1}{N_k^{(t)}} \sum_{i=1}^n r_{i,k}^{(t)} \boldsymbol{x}_i, \\
-    \boldsymbol{m}_{n,k}^{(t+1)} &= \frac{\kappa_0\boldsymbol{\mu}_0 + N_k^{(t)} \bar{\boldsymbol{x}}_k^{(t)}}{\kappa_0 + N_k^{(t)}}, \\
+    \boldsymbol{m}_{n,k}^{(t+1)} &= \frac{\kappa_0\boldsymbol{m}_0 + N_k^{(t)} \bar{\boldsymbol{x}}_k^{(t)}}{\kappa_0 + N_k^{(t)}}, \\
     \kappa_{n,k}^{(t+1)} &= \kappa_0 + N_k^{(t)}, \\
-    (\boldsymbol{W}_{n,k}^{(t+1)})^{-1} &= \boldsymbol{W}_0^{-1} + \sum_{i=1}^{n} r_{i,k}^{(t)} (\boldsymbol{x}_i-\bar{\boldsymbol{x}}_k^{(t)})(\boldsymbol{x}_i-\bar{\boldsymbol{x}}_k^{(t)})^\top + \frac{\kappa_0 N_k^{(t)}}{\kappa_0 + N_k^{(t)}}(\bar{\boldsymbol{x}}_k^{(t)}-\boldsymbol{\mu}_0)(\bar{\boldsymbol{x}}_k^{(t)}-\boldsymbol{\mu}_0)^\top, \\
+    (\boldsymbol{W}_{n,k}^{(t+1)})^{-1} &= \boldsymbol{W}_0^{-1} + \sum_{i=1}^{n} r_{i,k}^{(t)} (\boldsymbol{x}_i-\bar{\boldsymbol{x}}_k^{(t)})(\boldsymbol{x}_i-\bar{\boldsymbol{x}}_k^{(t)})^\top + \frac{\kappa_0 N_k^{(t)}}{\kappa_0 + N_k^{(t)}}(\bar{\boldsymbol{x}}_k^{(t)}-\boldsymbol{m}_0)(\bar{\boldsymbol{x}}_k^{(t)}-\boldsymbol{m}_0)^\top, \\
     \nu_{n,k}^{(t+1)} &= \nu_0 + N_k^{(t)},\\
     \alpha_{n,k}^{(t+1)} &= \alpha_{0,k} + N_k^{(t)}, \\
     \ln \rho_{i,k}^{(t+1)} &= \psi (\alpha_{n,k}^{(t+1)}) - \psi ( {\textstyle \sum_{k=1}^K \alpha_{n,k}^{(t+1)}} ) \notag \\

bayesml/gaussianmixture/test.py

@@ -1,5 +1,8 @@
 from bayesml import gaussianmixture
 import numpy as np
 
-model = gaussianmixture.LearnModel(num_classes=3, degree=2, h0_w_mats=np.identity(2)*2)
-print(model.calc_vl())
+model = gaussianmixture.LearnModel(num_classes=5, degree=3)
+
+x = np.random.rand(10,3)
+
+model.update_posterior(x,num_init=3)

doc/devdoc/__init__関数の仕様変更.md

@@ -0,0 +1,31 @@
+# __init__関数の仕様変更についての相談とお願い
+
+## 背景
+
+* ベイズ統計学的には事前分布と事後分布のハイパーパラメータは同じ形の方がうれしい．
+  * 事前分布の意味
+  * 逐次更新
+  * 混合正規分布の場合
+* 事前分布のハイパーパラメータを一般化すると，インスタンス生成時のハイパーパラメータ指定方法が多様になりすぎる．
+  * 個々の混合要素のハイパーパラメータを個別に指定する
+  * 全混合要素のハイパーパラメータを共通にする（ブロードキャストできると嬉しい）
+  * 次元がそろわない入力を受け付けなければならない．
+* __init__の実装の手間が大きい．
+
+## 対応案
+
+モデルのパラメータを3つに分ける．
+
+* constants
+  * 事後分布更新などに関与しないが，paramsやh_paramsの行列サイズを定める既知定数．degree, num_classesなど．
+  * インスタンス生成時，必ず手動で与える形式にする（デフォルト値無し）
+  * c_valname
+* params
+  * サンプルを生成する確率分布のパラメータ
+* h_params, h0_params, hn_params
+  * paramsを生成する確率分布（事前分布）のパラメータ
+  * これらは同じ型の行列にする
+* p_params
+  * 予測分布のパラメータ
+
+次元確認は利用者が手動で与えたconstantsに対してparamsやh_paramsが整合しているかのみをチェックする．その際，ブロードキャストもうまく利用する．