LoginSignup
11
12

More than 5 years have passed since last update.

@yh0sh

ノイズ有り対応トピックモデルの実装

Last updated at Posted at 2016-07-23

本記事を読むための事前知識

本記事では通常のトピックモデルについての説明は省きます.
通常のトピックモデルについては,Qiitaでは以下の記事が理解しやすいと思います.
数式をなるべく使わずにトピックモデルの解説にチャレンジ

実装まで挑戦したいという方には,以下のブログがとても参考になります(私も参考にさせて頂きました).
Latent Dirichlet Allocations の Python 実装 (Mi manca qualche giovedi`?)

ノイズ有り対応トピックモデル

トピックモデルの拡張として,補助情報についてもトピックを区分する対応トピックモデルが2003年に発表されました.
補助情報とは,記事のタグや商品のカテゴリーなど,文書の内容に対して補助的に与えられる情報のことです.
対応トピックモデルでは,全ての補助情報がトピックに割り当てられますが,文書の内容と関係しない一般補助情報も同様にトピックが割り当てられます.一般補助情報とは,例えばSNSでは「いいね」や「後で読む」,Amazonでは「☆☆☆☆☆」(評価点)といった記事や商品の内容とは関係しない語彙のことです.一般補助情報をノイズとして,トピックとは別に区分するトピックモデルをノイズ有り対応トピックモデルといいます1.

生成過程

ノイズ有り対応トピックモデルの生成過程では，通常のトピックモデルに補助情報分布，トピック毎の補助情報の生成過程が加えています．
トピック毎の補助情報分布は，文書$d$内の単語$m$を生成したトピック$k$の割合($N_{dk}/N_{d}$)に応じて割り当てられ，また，補助情報と一般補助情報を区別する二値の潜在変数($r_{dm} ∈\{0,1\}$)によって，補助情報が補助情報分布から生成されるのか($r_{dm} = 1$のとき)，一般補助情報分布から生成されるのか($r_{dm} = 0$のとき)を決定することが，ノイズ有り対応トピックモデルの特徴です．

以下の観測変数$\{(w_{d},t_{d})\}^{D}_{d=1}$の生成過程によって,モデルを生成します.

1. 関係性分布を生成 $λ$ ~ $Beta(η)$

2. 一般補助情報分布を生成 $Ψ_{0}$~$Dirichlet(γ)$

3. For トピック$k=1,...,K$

(a) 単語分布を生成 $φ_{k}$~$Dirichlet(β)$

(b) 補助情報分布を生成 $Ψ_{k}$~$Dirichlet(γ)$

4. For 文書$d=1,...,D$

(a) 単語トピック分布を生成 $Θ_{d}$~$Dirichlet(α)$

(b) For 単語数 $n=1,...,N_{d}$

i. トピック$z_{dn}$~$Categorical(Θ_{d})$

ii. 単語$w_{dn}$~$Categorical(φ_{zdn})$

(c) For 補助情報数 $m=1,...,M_{d}$

i. トピック $c_{dm}$~$Categorical( \{ \frac{N_{kd}}{N_{d}} \}^{K}_{k=1})$

ii. 関係性 $r_{dm}$~$Bernoulli(λ)$

iii. 補助情報 $t_{dm}$

$\space Categorical(Ψ_{0})\space$if $rdm = 0$

$\space Categorical(Ψ_{cdm})\space$otherwise

グラフィカルモデル表現

ノイズ有り対応トピックモデルのグラフィカルモデル表現は以下の図で示されます.補助情報ではトピック数(K) + 1 回の繰り返しによって,補助情報のトピックと一般補助情報を区分します.

graphical_model.jpg

記号 説明
$w$ 単語(観測変数)
$t$ 補助情報(観測変数)
$z$ 単語トピック
$θ$ 単語トピック分布
$Φ$ 単語分布
$c$ 補助情報トピック
$ψ$ 補助情報分布
$r$ 関係性
$λ$ 関係性分布
$α$ $θ$を生成するディリクレ分布のパラメータ
$β$ $Φ$を生成するディリクレ分布のパラメータ
$γ$ $ψ$を生成するディリクレ分布のパラメータ
$η$ $λ$を生成するベータ分布のパラメータ
$D$ 文書数
$N$ ある文書の中で現れる単語数
$M$ ある文書の中で現れる補助情報数
$K$ トピック数
$V$ 全文書の中で現れる単語の種類数(語彙数)

サンプリング確率

グラフィカルモデル表現で示される，未知変数$z$，$c$，$r$の集合$Z$，$C$，$R=(r_{11},...,r_{DN_{d}})$のサンプリング確率は，崩壊型ギプスサンプリングによって以下の式で示されます($r_{DN_{d}}$は文書$D$にある$N_{d}$番目の単語の関係性)．

p(z_{dn}=k|W,X,Z_{\dn},C,R,α,β,γ,η)
∝(N_{dk\dn}+α)\frac{N_{kw_{dn}\dn}+β}{N_{k\dn}+βV} \left( \frac{N_{dk\dn}+1}{N_{dk\dn}} \right)^{M_{dk}}

p(c_{dm}=k|W,X,Z,C_{\dm},R,α,β,γ,η)
∝N_{dk}\frac{M_{kx_{dm}\dm}+γ} {M_{k\dm}+γS}

p(r_{dm}=0|W,X,Z,C,R_{\dm},α,β,γ,η)
∝(M_{0\dm}+η)\frac{M_{0x_{dm}\dm}+γ} {M_{0\dm}+γS}

p(r_{dm}=1|W,X,Z,C,R_{\dm},α,β,γ,η)
∝(M_{1\dm}+η)\frac{M_{c_{dm}x_{dm}\dm}+γ} {M_{c_{dm}\dm}+γS}
記号 説明
$N_{dk}$ 文書$d$でトピック$k$が割り当てられた単語数
$N_{kv}$ 文書集合全体で語彙$v$にトピック$k$が割り当てられた単語数
$N_{k}$ 文書集合全体でトピック$k$が割り当てられた単語数
$M_{0}$($M_{1}$) 内容と関係ない(ある)補助情報の数
$M_{0s}$ 内容と関係ない補助情報$s$の数
$M_{dk}$ 文書$d$でトピック$k$が割り当てられた補助情報の数
$M_{ks}$ 内容と関係があり，かつトピック$k$が割り当てられた補助情報$s$の数
$M_{k}$ 文書集合全体でトピック$k$が割り当てられた補助情報の数

※\dn：文書$d$の$n$番目の単語を除いた時の単語数を指します．例えば，$N_{dk\dn}$は$n$番目の単語を除いた時の文書$d$でトピック$k$が割り当てられた単語数となります．
※\dm：文書$d$の$m$番目の補助情報を除いた時の補助情報の数を指します．

データセットについて

今回，Qiitaの最新記事9999本(2016/7/17時点)から抽出した単語と，補助情報にタグを使用したデータセットを作成しました．

Qiitaの記事とタグ取得については，以下の記事を参考にさせて頂きました．
Qiitaの投稿をガーッと取得するバッチ処理用QiitaAPI Pythonラッパー

MeCabを使用した形態素解析では，以下の記事を参考にさせて頂きました．
PythonでWebスクレイピングと形態素解析。

形態素解析の結果より，1本の記事の中で頻出回数の多い上位30位の単語を使用します．

今回のデータセットはソースコードと共にGitHubにアップロードしています2 3．

実装

Python 2.7を使用して,ノイズ有り対応トピックモデルを実装しました.モジュールはNumpyのみで行えます.
解法アルゴリズムとしては,トピックモデルの解法ではお馴染みのGibbs Samplingを用います.
なお，Numpyのみでは，配列要素探索に時間がかかりすぎてしまうため，Cython+Numpyを用いてチューニングしています．

corr_lda.pyx
# distutils: language=c++
# -*- coding: utf-8 -*-

import os
import math
from scipy import special 
from d_profile import profile

import numpy as np
cimport numpy as np
from libcpp.vector cimport vector

DTYPE = np.int
FTYPE = np.float
ctypedef np.int_t DTYPE_t
ctypedef np.float_t FTYPE_t

cdef class CorrLDA:
  # MAIN VARIABLES
  cdef int K, D, V
  cdef double alpha, beta
  cdef np.ndarray Nk
  cdef np.ndarray Ndk
  cdef np.ndarray Nkv
  cdef np.ndarray Nd
  cdef np.ndarray zdn # topic sequence for word n in document d
  cdef np.ndarray vdn # vocaburary sequence for word n in document d
  # SUPPLY INFORMATION VARIABLES
  cdef int S
  cdef double M0, M1
  cdef double gamma, eta, lam
  cdef np.ndarray Mk
  cdef np.ndarray Mdk
  cdef np.ndarray M0s
  cdef np.ndarray Mks
  cdef np.ndarray Md
  cdef np.ndarray xdm # topic sequence for supply m in document d
  cdef np.ndarray sdm # vocaburary sequence for supply m in document d
  cdef np.ndarray rdm # relevance sequence for supply m in document d

  def __init__(self, K, alpha, beta, eta, lam):
    self.K = K
    self.alpha = alpha
    self.beta = beta
    self.eta = eta
    self.lam = lam
    self.gamma = np.random.beta(self.lam, self.lam)

  # make document dictionary with a key of document id
  def make_dict(self, vfile):
    vdict = {}
    fvoc = open(vfile)
    for num,line in enumerate(fvoc):
      vdict[num] = line[:-1]
    fvoc.close()
    return vdict

  # count up Nk, Ndk, Nkv, and Nd
  def set_corpus(self, dfile, prof = False):
    self.Nk = np.zeros(self.K,dtype=FTYPE)
    fdoc = open(dfile)
    pre_doc = 0
    cdef vector[int] zdn
    cdef vector[int] vdn
    cdef int count
    for num, line in enumerate(fdoc):
      if num == 0:
        self.D = int(line[:-1])
        self.Ndk = np.zeros([self.D, self.K],dtype=FTYPE)
        self.Nd = np.zeros(self.D,dtype=DTYPE)
        continue
      elif num == 1:
        self.V = int(line[:-1])
        self.Nkv = np.zeros([self.K, self.V],dtype=FTYPE)
        continue
      elif num == 2:
        continue
      else:
        d, v, c = line[:-1].split(' ')
        k = np.random.randint(self.K,size=1)[0]
        doc = int(d) - 1
        voc = int(v) - 1
        count = int(c)
        self.Ndk[doc, k] += count
        self.Nkv[k, voc] += count
        self.Nk[k] += count
        self.Nd[doc] += count
        for c in range(count):
          vdn.push_back(voc)
          zdn.push_back(k)
        if pre_doc != doc and prof:
          print pre_doc+1, self.Ndk[pre_doc], self.Nd[pre_doc]
          pre_doc = doc
    if prof:
      print pre_doc+1, self.Ndk[pre_doc], self.Nd[pre_doc]
    fdoc.close()
    self.Nk += self.beta * self.V
    #self.Ndk += self.alpha
    self.Nkv += self.beta
    self.zdn = np.array(list(zdn))
    self.vdn = np.array(list(vdn))
    return 0

  # count up M0, M1, Mk, Mdk, M0s, Mks, and Md
  def set_supply(self, dfile, prof = False):
    self.M0 = 0
    self.M1 = 0
    self.Mk = np.zeros(self.K + 1,dtype=FTYPE) # K + 1 includes the index of noisy topic
    self.Mdk = np.zeros([self.D, self.K],dtype=FTYPE)
    self.Md = np.zeros(self.D,dtype=DTYPE) # for numbering iteration only, not using calculation
    fdoc = open(dfile)
    pre_doc = 0
    cdef vector[int] xdm
    cdef vector[int] sdm
    cdef vector[int] rdm
    cdef int count, num_noisy
    for num, line in enumerate(fdoc):
      if num == 0 or num == 2:
        continue
      elif num == 1:
        self.S = int(line[:-1])
        self.Mks = np.zeros([self.K, self.S],dtype=FTYPE)
        self.M0s = np.zeros(self.S,dtype=FTYPE)
        continue
      else:
        r = np.random.binomial(1, self.gamma) # setting relevance randomly
        d, s, c = line[:-1].split(' ')
        doc = int(d) - 1
        supply = int(s) - 1
        count = int(c)
        if r == 1: # relevance
          k = np.random.multinomial(1,self.Ndk[doc] / self.Ndk[doc].sum()).argmax()
          self.Mdk[doc, k] += count
          self.Mks[k, supply] += count
          self.M1 += count
        else: # irrelevance
          k = self.K  
          self.M0s[supply] += count
          self.M0 += count
          num_noisy += count
        self.Mk[k] += count
        self.Md[doc] += count
        for c in range(count):
          sdm.push_back(supply)
          xdm.push_back(k)
          rdm.push_back(r)
        # if prof is true, print out below:
        # doc-id, word count of all topics(array), that of noisy topic, total count of words in the document with doc-id
        if pre_doc != doc and prof:
          print pre_doc+1, self.Mdk[pre_doc], num_noisy, self.Md[pre_doc]
          pre_doc = doc
          num_noisy = 0
    if prof:
      print pre_doc+1, self.Mdk[pre_doc], num_noisy, self.Md[pre_doc]
    fdoc.close()
    self.M0 += self.eta
    self.M1 += self.eta
    self.Mk += self.gamma * self.S
    self.Mks += self.gamma
    self.M0s += self.gamma
    self.xdm = np.array(list(xdm))
    self.sdm = np.array(list(sdm))
    self.rdm = np.array(list(rdm))
    return 0

  # if prof is true, print out the profile of the inference
  @profile
  def inference_prof(self):
    return self.inference()

  def inference(self):
    # MAIN VARIABLES
    cdef int t, v, new_t
    cdef int idx = 0 # index of topic for main
    cdef np.ndarray[FTYPE_t, ndim=1] theta = np.zeros(self.K,dtype=FTYPE)
    cdef np.ndarray[FTYPE_t, ndim=1] phi = np.zeros(self.K,dtype=FTYPE)
    cdef np.ndarray[FTYPE_t, ndim=1] p_z = np.zeros(self.K,dtype=FTYPE)
    # SUPPLY INFORMATION VARIABLES
    cdef int k, s, new_ts, ktemp = 0
    cdef int idxs = 0 # index of topic for supply information
    cdef double p_r0
    cdef np.ndarray[FTYPE_t, ndim=1] p_r1 = np.zeros(self.K,dtype=FTYPE)
    cdef np.ndarray[FTYPE_t, ndim=1] p_y = np.zeros(self.K,dtype=FTYPE)

    for d in range(len(self.Nd)):
      # MAIN ITERATION
      for n in range(self.Nd[d]):
        t = self.zdn[idx]
        v = self.vdn[idx]
        # REDUCE COUNT OF OLD TOPIC
        self.Ndk[d][t] -= 1
        self.Nkv[t][v] -= 1
        self.Nk[t] -= 1
        # INFERENCE
        theta = self.Ndk[d] + self.alpha
        phi = self.Nkv[:,v] / self.Nk
        p_z = theta * phi * ((self.Ndk[d] + 1.0) / self.Ndk[d])**self.Mdk[d]
        new_t = np.random.multinomial(1,p_z / p_z.sum()).argmax() # decide new topic as categorical dist
        # COUNT ADD TO NEW TOPIC
        self.zdn[idx] = new_t
        self.Ndk[d][new_t] += 1
        self.Nkv[new_t][v] += 1
        self.Nk[new_t] += 1
        idx += 1

      # ITERATION FOR SUPPLY INFOMATION
      for n_s in range(self.Md[d]):
        k = self.xdm[idxs] # old topic
        s = self.sdm[idxs]
        # REDUCE COUNT OF OLD TOPIC
        self.Mk[k] -= 1
        if self.rdm[idxs] == 1:
          self.Mdk[d, k] -= 1
          self.Mks[k, s] -= 1
          self.M1 -= 1
        else:
          self.M0s[s] -= 1
          self.M0 -= 1
        p_r1 = self.M1 * self.Mks[:, s] / self.Mk[:self.K]
        p_r0 = self.M0 * self.M0s[s] / self.Mk[self.K]
        r = np.random.binomial(1, p_r1.sum() / (p_r0 + p_r1.sum()))
        # COUNT ADD TO NEW TOPIC
        self.rdm[idxs] = r # replace to new rdm = {0,1}
        if r == 1:
          p_y = self.Ndk[d] * (self.Mks[:,s] + self.gamma) / self.Mk[k]
          new_ts = np.random.multinomial(1,p_y / p_y.sum()).argmax()
          self.xdm[idxs] = new_ts # replace to new topic
          self.Mdk[d, new_ts] += 1
          self.Mks[new_ts, s] += 1
          self.M1 += 1
        else:
          self.xdm[idxs] = self.K # replace to noisy topic
          self.M0s[s] += 1
          self.M0 += 1
        self.Mk[self.xdm[idxs]] += 1
        idxs += 1
    cdef double perplexity = self.__calc_perplexity()
    return perplexity

  # calculation perplexity per one iteration
  def __calc_perplexity(self):
    cdef int s, N, idx = 0
    cdef double log_per = 0.0
    cdef double Kalpha = self.K * self.alpha
    cdef np.ndarray[FTYPE_t, ndim=1] theta = np.zeros(self.K,dtype=FTYPE)
    cdef np.ndarray[FTYPE_t, ndim=2] phi = self.Nkv / self.Nk[:,np.newaxis]
    for d in range(len(self.Nd)):
      theta = self.Ndk[d,:] / (self.Nd[d] + Kalpha) * ((self.Ndk[d] + 1.0) / self.Ndk[d])**self.Mdk[d]
      for v in range(self.Nd[d]):
        v = self.vdn[idx]
        log_per -= np.log(np.inner(phi[:,v], theta))
        idx += 1
    N = idx + 1
    return np.exp(log_per / N)

結果

トピック数は50として，反復計算200回行ったときの結果を示します．

トピック例

トピック1-5までの単語とタグ，そしてノイズとなったタグを下表に示します．
各トピックの中で確率の高い上位10位の単語とタグ(補助情報)，それらのトピック内の割合を順に並べています．

トピック毎に関連単語，関連タグが集まっていますが，同一トピック内の単語と補助情報の関連性については微妙です．
タグの種類が少ないため，ノイズトピックも含めて，複数のトピックに同じタグが存在しています．
ある一つのトピックだけに関連するタグは，あまり見られないということだと思います．

topic 単語 タグ
topic 1 com : 0.0143840966187 iOS : 0.0137813506742
よう : 0.0118339832479 JavaScript : 0.0104194960757
https : 0.0110616631985 Ruby : 0.0102093801633
http : 0.00810353168844 Windows : 0.00915880060132
こと : 0.00746236032665 Mac : 0.00873856877652
場合 : 0.00708348634014 Android : 0.00852845286412
の : 0.00673375650643 Rails : 0.00810822103931
以下 : 0.00552427416487 Linux : 0.0072677573897
設定 : 0.00542226963004 Swift : 0.00663740965249
タスク : 0.00447508466376 Python : 0.00642729374009
topic 2 tf : 0.038923153859 Python : 0.0185783855201
train : 0.0312198497729 TensorFlow : 0.0176585006859
test : 0.0220516550738 DeepLearning : 0.0136723330711
data : 0.0185914824188 Chainer : 0.00692651095376
model : 0.0179095505817 機械学習 : 0.00447348472926
loss : 0.0171518485404 MachineLearning : 0.00447348472926
np : 0.0124793526194 Keras : 0.0041668564512
step : 0.0115953669046 深層学習 : 0.0041668564512
accuracy : 0.00944854445438 CNN : 0.00386022817314
print : 0.00901917996434 ConvolutionalNeuralNetworks : 0.00355359989508
topic 3 mysql : 0.0438792621576 MySQL : 0.0172747954451
el : 0.0398385043809 VBA : 0.0126760035933
As : 0.0384579121405 Excel : 0.00479236041866
的 : 0.0315886239199 ExcelVBA : 0.00282144962501
中 : 0.0257968711066 access : 0.00216447936046
Long : 0.0166378201459 Xcode : 0.00150750909591
MySQL : 0.0123613514988 VB.Net : 0.00150750909591
remi : 0.0118225837952 Azure : 0.00117902396363
Sub : 0.0117552378323 mecab : 0.00117902396363
Function : 0.0114858539805 mariadb : 0.00117902396363
topic 4 stocks : 0.133028236329 WebAPI : 0.000907792530452
jp : 0.10110874249 努力1mm : 0.000907792530452
http : 0.0886481914183 #slick-codegen : 0.000557195545988
日 : 0.08622056503 openssl : 0.000557195545988
code : 0.0815739363962 AdobeAIR : 0.000557195545988
co : 0.0778755993203 mex : 0.000557195545988
円 : 0.0719013625053 日時 : 0.000557195545988
detail : 0.0682409570917 zybo : 0.000557195545988
yahoo : 0.0671409388845 Carbon : 0.000557195545988
株 : 0.0664961006251 Chrome : 0.000557195545988
topic 5 data : 0.115147263212 crystal : 0.00490350058429
value : 0.0738510543501 Ruby : 0.00322298524599
key : 0.0632111293615 riot : 0.00120636684003
sample : 0.0343835825957 ぐるナビ : 0.00120636684003
fa : 0.0293961177574 Ansible : 0.00120636684003
default : 0.0268691355726 Qiita : 0.00120636684003
size : 0.0263371393232 PHP : 0.00120636684003
class : 0.0221144190933 Twitter : 0.00120636684003
option : 0.0123389880101 websocket : 0.000870263772375
chef : 0.0120729898854 Parslet : 0.000870263772375
noisy topic Go : 0.000211416490486
TensorFlow : 0.000211416490486
Vim : 0.000211416490486
C# : 0.000211416490486
めざせライブラリマスター : 0.000211416490486
Ubuntu : 0.000211416490486
機械学習 : 0.000211416490486
Windows : 0.000211416490486
PostgreSQL : 0.000211416490486
VirtualBox : 0.000211416490486

Perplexity

毎回の反復計算でトピックモデルのperplexityを求めています．
50回前後でperplexityが一定となった確率モデルが生成されたことを示しています．

perplexity2.png

ソースコード

全てのソースコード，サンプルデータはgithubにアップロードしています．

参考文献

  1. 正確には,ノイズ有り対応トピックモデルは,対応トピックモデルの拡張です. ↩

  2. データセットのファイル構成はUCI Machine Learning Repository: Bag of Words Data Setに準じています． ↩

  3. 驚くほど拙い英文のREADMEが付いていますので失笑しながらご覧ください．．． ↩

11
12
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
11
12