唯物是真 @Scaled_Wurm

プログラミング(主にPython2.7)とか機械学習とか

Twitterの投稿時間で類似度を計算してみた - 確率分布の類似度

以前集合やベクトルの類似度の記事を書いたんですが、確率分布の類似度には触れていなかったのでついでに書きました

ツイート時間分布の類似度を求める

今回はツイート時間ごとの頻度を正規化して、確率分布とみなして類似度を計算してみます
私のアカウント(以下mainと表記)に対して、私のもう一つのアカウント(以下subと表記)+私がリプライを送ってる数が多い(以下friendと表記)上位5人と比較します
subがfriendよりも似た結果になることを期待しています

以下にツイート時間の分布(main + sub + friend 5人)を載せました
ある程度似ていますが、人によってそこそこ形が違っていて、特に午前中の投稿時間の差は特徴的に見えます
またfriend1の一人だけは大きく違った傾向を示しています
f:id:sucrose:20131107121242p:plain
以下ではこれらが定量的にどれぐらい異なるのか類似度を計算して比較していきます

グラフを描くのに使ったPythonソースコード(これ以降掲載するコードはnumpyやscipyなどのライブラリを使っています

import numpy as np
from pylab import *

data = []
#main
data.append(np.array([4006,2064,856,385,162,112,163,234,360,871,1657,2435,2903,2626,2790,2406,2366,2523,2841,2552,2803,2464,3421,4457], dtype='float'))
#sub
data.append(np.array([593,418,199,65,23,8,7,28,31,90,236,317,369,403,262,231,186,178,220,305,340,368,547,633], dtype='float'))
#friend 1
data.append(np.array([20,10,9,5,9,12,36,338,121,42,32,34,162,98,35,32,54,55,65,203,179,159,77,38], dtype='float'))
#friend 2
data.append(np.array([7458,3890,1293,323,235,280,673,2169,4569,5515,4976,4513,4225,4150,3416,3047,3179,3642,4296,4911,4278,5158,7282,7939], dtype='float'))
#friend 3
data.append(np.array([1254,803,307,176,48,61,231,512,528,554,750,704,702,549,575,533,560,549,456,553,759,922,1033,1317], dtype='float'))
#friend 4
data.append(np.array([816,382,161,23,8,23,61,227,460,621,523,603,820,730,596,580,569,619,596,770,792,857,950,937], dtype='float'))
#friend 5
data.append(np.array([2219,1429,912,479,256,189,246,370,624,691,921,878,1362,1053,879,950,1187,1039,1109,1221,1495,1580,2237,2236], dtype='float'))
#normalize
data_prob = [d / sum(d) for d in data]

xlim(0, 23)
plot(np.array(data_prob)[:7,:].T, linewidth=5)
legend('main sub friend1 friend2 friend3 friend4 friend5'.split())

show()

使ったデータ

過去のツイートを記録してくれるWebサービスであるTwilog (私のアカウント)では、アカウントごとの様々な統計情報を表示してくれます
f:id:sucrose:20131107015953p:plain

それらの統計情報の一部はグラフで表示されていますが、実は以下のようにtxtで数値が簡単に入手できたりします
ただし統計のページにアクセスした直後にしか存在しなさそう(?)

ディレクトリ名の"s"の部分はおそらくアカウント名の頭文字になっていて、末尾の"4"などの部分がデータの種類になっています
このtxtのファイル名はグラフの表示に使われているFlashに与えられているのでブラウザのツールなどでファイル名を見れば簡単にわかります

確率分布の類似度

簡単に確率分布の類似度について説明します
記事の冒頭から類似度といってきましたが、確率分布についてはどちらかというと違いの大きさを表している尺度(距離のようなもの)といったほうが正しいです

特によく使われているKLダイバージェンスと、その仲間であるJSダイバージェンスを紹介します
以下の数式では対数の底を2としていますが、自然対数などの他の底を使っても単位や値域が変わるだけでだいたい一緒です

KLダイバージェンス(Kullback–Leibler divergence)

KLダイバージェンス情報理論に基づいた確率分布の違いの大きさを表す尺度で、エントロピー(情報量)と関係があります
他にも情報利得(Information gain)とか様々な呼び名で呼ばれます

KLダイバージェンスの式は離散確率分布\(P, Q\)について次のように定義されています
$$D_{\mathrm{KL}}(P\|Q) = \sum_i P(i) \log_2 \frac{P(i)}{Q(i)}$$
ただし\(0 \log_2 0 = 0\)とします

\(D_{\mathrm{KL}}(P\|Q)\)は\(0\)以上の値をとり、大きければ大きいほど分布の違いが大きいことを表しています
注意点として、KLダイバージェンスには対称性がない(\(D_{\mathrm{KL}}(P\|Q) \neq D_{\mathrm{KL}}(Q\|P)\))こと、\(Q(i)=0\)の場合には有限の値にならないことがあります

計算自体は単純で簡単に計算できます
例えばコインの表裏(H, T)が出る確率がそれぞれ\([\frac{1}{2}, \frac{1}{2}]\)の確率分布を\(P\)、\([\frac{2}{3}, \frac{1}{3}]\)の確率分布を\(Q\)とおきます
すると、KLダイバージェンスは以下のようになります
$$
\begin{eqnarray}
D_{\mathrm{KL}}(P\|Q) &=& \sum_i P(i) \log_2 \frac{P(i)}{Q(i)} = P(H) \log_2 \frac{P(H)}{Q(H)} + P(T) \log_2 \frac{P(T)}{Q(T)}\\
&=& \frac{1}{2} \log_2 \frac{\frac{1}{2}}{\frac{2}{3}} + \frac{1}{2}\log_2 \frac{\frac{1}{2}}{\frac{1}{3}} = \frac{1}{2} \log_2 \frac{3}{4} + \frac{1}{2}\log_2 \frac{3}{2}\\
&\sim& 0.085
\end{eqnarray}
$$
関数もnumpyを使えば1行で書けます(リスト内包表記とかを使えばnumpyなしでも書けそうですが

import numpy as np
def kl(p, q):
    return np.sum(np.where(p != 0, p * np.log2(p / q), 0))

過去のいくつかの記事にも登場しているように、この尺度は機械学習を勉強しているとよくでてきます
PRMLなどを読むと最尤推定などのモデルのパラメータの最適化が結局KLダイバージェンスの最小化と等しい、という話になることが多いです

パターン認識と機械学習 上

パターン認識と機械学習 上

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

2つの確率分布\(P\)と\(Q\)があったとき、\(P\)と\(Q\)のKLダイバージェンスは\(P\)と\(Q\)のクロスエントロピーと\(P\)のエントロピーの差に等しくなります
ちなみにクロスエントロピーも確率分布の違いを表す尺度として使われていて、言語モデルが実際のテキストに対してどれぐらいうまくいっているのかをはかるのに使われているのをみかけます

データで示した各ユーザーのツイート時間分布を使って、mainに対するsub + friend 5人との間のKLダイバージェンスを求めました

明らかに形の違うfriend 1がもっとも大きく、またsubが一番小さなKLダイバージェンスとなっており期待通りの結果になっています
f:id:sucrose:20131108022502p:plain

import numpy as np
from pylab import *

def kl(p, q):
    return np.sum(np.where(p != 0, p * np.log2(p / q), 0))

data = []
#main
data.append(np.array([4006,2064,856,385,162,112,163,234,360,871,1657,2435,2903,2626,2790,2406,2366,2523,2841,2552,2803,2464,3421,4457], dtype='float'))
#sub
data.append(np.array([593,418,199,65,23,8,7,28,31,90,236,317,369,403,262,231,186,178,220,305,340,368,547,633], dtype='float'))
#friend 1
data.append(np.array([20,10,9,5,9,12,36,338,121,42,32,34,162,98,35,32,54,55,65,203,179,159,77,38], dtype='float'))
#friend 2
data.append(np.array([7458,3890,1293,323,235,280,673,2169,4569,5515,4976,4513,4225,4150,3416,3047,3179,3642,4296,4911,4278,5158,7282,7939], dtype='float'))
#friend 3
data.append(np.array([1254,803,307,176,48,61,231,512,528,554,750,704,702,549,575,533,560,549,456,553,759,922,1033,1317], dtype='float'))
#friend 4
data.append(np.array([816,382,161,23,8,23,61,227,460,621,523,603,820,730,596,580,569,619,596,770,792,857,950,937], dtype='float'))
#friend 5
data.append(np.array([2219,1429,912,479,256,189,246,370,624,691,921,878,1362,1053,879,950,1187,1039,1109,1221,1495,1580,2237,2236], dtype='float'))
#normalize
data_prob = [d / sum(d) for d in data]

result = [kl(data_prob[0], data_prob[i]) for i in xrange(1, len(data_prob))]
bar(xrange(len(result)), result, 0.8, align='center', color='grcmyk')
xticks([i for i in xrange(len(result))], 'sub friend1 friend2 friend3 friend4 friend5'.split())
show()

JSダイバージェンス(Jensen–Shannon divergence)

JSダイバージェンスは以下のように定義されています
$$D_{\mathrm{JS}} = \frac{1}{2} D_{\mathrm{KL}} \left (P \| M \right ) + \frac{1}{2} D_{\mathrm{KL}}\left (Q \| M \right )$$
ただし\(M = \frac{1}{2}(P+Q)\)

\(M\)の式を見ればわかるように、\(P\)と\(Q\)の平均的な分布との間のKLダイバージェンスを求めています
またKLダイバージェンスとは異なり対称性があり、\(Q(i)\)が\(0\)でも有限の値になるという扱いやすい性質があります
JSダイバージェンス平方根距離の定義を満たしているそうです
値域は底が2の対数を使った時には\(0 \leq D_{\mathrm{JS}} \leq 1\)

以下のようにKLダイバージェンスと似た結果が得られました
f:id:sucrose:20131108022952p:plain

import numpy as np
from pylab import *

def kl(p, q):
    return np.sum(np.where(p != 0, p * np.log2(p / q), 0))

def js(p, q):
    m = (p + q) / 2
    return (kl(p, m) + kl(q, m)) / 2.0

data = []
#main
data.append(np.array([4006,2064,856,385,162,112,163,234,360,871,1657,2435,2903,2626,2790,2406,2366,2523,2841,2552,2803,2464,3421,4457], dtype='float'))
#sub
data.append(np.array([593,418,199,65,23,8,7,28,31,90,236,317,369,403,262,231,186,178,220,305,340,368,547,633], dtype='float'))
#friend 1
data.append(np.array([20,10,9,5,9,12,36,338,121,42,32,34,162,98,35,32,54,55,65,203,179,159,77,38], dtype='float'))
#friend 2
data.append(np.array([7458,3890,1293,323,235,280,673,2169,4569,5515,4976,4513,4225,4150,3416,3047,3179,3642,4296,4911,4278,5158,7282,7939], dtype='float'))
#friend 3
data.append(np.array([1254,803,307,176,48,61,231,512,528,554,750,704,702,549,575,533,560,549,456,553,759,922,1033,1317], dtype='float'))
#friend 4
data.append(np.array([816,382,161,23,8,23,61,227,460,621,523,603,820,730,596,580,569,619,596,770,792,857,950,937], dtype='float'))
#friend 5
data.append(np.array([2219,1429,912,479,256,189,246,370,624,691,921,878,1362,1053,879,950,1187,1039,1109,1221,1495,1580,2237,2236], dtype='float'))
#normalize
data_prob = [d / sum(d) for d in data]

result = [js(data_prob[0], data_prob[i]) for i in xrange(1, len(data_prob))]
bar(xrange(len(result)), result, 0.8, align='center', color='grcmyk')
xticks([i for i in xrange(len(result))], 'sub friend1 friend2 friend3 friend4 friend5'.split())
show()

その他

他にも様々な尺度がありますが紹介するのは上の2つにとどめておくので以下のWikipediaの記事などを参考にしてください

今回は確率分布の類似度としてダイバージェンスを計算しましたが、単純にベクトルとみなしてコサイン類似度などを求めてもある程度似た結果が得られると思います
以下コサイン類似度の結果
f:id:sucrose:20131108023311p:plain

import numpy as np
import scipy.spatial.distance
from pylab import *

data = []
#main
data.append(np.array([4006,2064,856,385,162,112,163,234,360,871,1657,2435,2903,2626,2790,2406,2366,2523,2841,2552,2803,2464,3421,4457], dtype='float'))
#sub
data.append(np.array([593,418,199,65,23,8,7,28,31,90,236,317,369,403,262,231,186,178,220,305,340,368,547,633], dtype='float'))
#friend 1
data.append(np.array([20,10,9,5,9,12,36,338,121,42,32,34,162,98,35,32,54,55,65,203,179,159,77,38], dtype='float'))
#friend 2
data.append(np.array([7458,3890,1293,323,235,280,673,2169,4569,5515,4976,4513,4225,4150,3416,3047,3179,3642,4296,4911,4278,5158,7282,7939], dtype='float'))
#friend 3
data.append(np.array([1254,803,307,176,48,61,231,512,528,554,750,704,702,549,575,533,560,549,456,553,759,922,1033,1317], dtype='float'))
#friend 4
data.append(np.array([816,382,161,23,8,23,61,227,460,621,523,603,820,730,596,580,569,619,596,770,792,857,950,937], dtype='float'))
#friend 5
data.append(np.array([2219,1429,912,479,256,189,246,370,624,691,921,878,1362,1053,879,950,1187,1039,1109,1221,1495,1580,2237,2236], dtype='float'))
#normalize
data_prob = [d / sum(d) for d in data]

result = [scipy.spatial.distance.cosine(data_prob[0], data_prob[i]) for i in xrange(1, len(data_prob))]
bar(xrange(len(result)), result, 0.8, align='center', color='grcmyk')
xticks([i for i in xrange(len(result))], 'sub friend1 friend2 friend3 friend4 friend5'.split())
show()

もしかしたら頻度に対する統計的検定とみなしてp値の比較を行ってもいいかもと考えて、試しにmainの分布に対するカイ二乗検定で適合度(goodness of fit)を調べてみました
しかし、p値が非常に小さくなってしまって意味のありそうな結果はまったく得られませんでした(以下の出力の右側がp値
そもそもp値同士を比較するのは検定として正しくなさそうな使い方の気がしますし、頻度の大小の影響も大きそうなのでそもそも向いてないのかもしれません(あまり検定には詳しくない

>>> scipy.stats.chisquare(data[1], data_prob[0] * sum(data[1]))
(475.68282464776814, 3.9976166273789829e-86)

まとめ

ツイート時間の確率分布を求めて類似度を計算した

Twilogからは統計情報が意外と簡単に入手できる(ユーザーが登録していれば

確率分布の違いの大きさはKLダイバージェンスなどを使うことで定量的にはかれる

もしも隠れて複垢をするときは、気をつけていないとツイート時間の分布が似ているかもしれないので、ツイート時間をずらすように気をつけましょう(そこまで見る人いないだろうけど

余談ですが\(0 \log 0\)の話に絡んで読んだ、\(0^0\)の記事が面白かったので紹介しておきます