隐含狄利克雷分布(Latent Dirichlet Allocation，LDA)
和机器学习中的线性判别分析(Linear Discriminant Analysis)不一样
首先是对狄利克雷分布的大致介绍
其就是对多项分布的先验分布.先验分布就是指对于可能出现的分布给出一个预设分布(根据背景知识),然后根据实际的样本信息去更新这个分布,最终得到后验分布.
beta分布是二项分布的先验分布,dirichlet分布则是多项分布的先验分布
从另一个方面来讲,狄利克雷分布就是分布的分布.也就是说,在使用多项分布时你所选用的多项分布也是由另一个分布选取的.
使用 AnkurTomar文章中的例子就是你有一批骰子,每一个骰子都是不均匀的,因此在进行掷骰子的实验时,骰子本身就是一个多项分布,选择哪一个骰子也是基于一个分布.这样最终得到的结果就是基于分布的分布也就是dirichlet分布.
LDA主题模型的一些假设前提
每一篇文章由多个主题的分布组成
每一个主题会包含一组词的分布
文章的每个词来自于一个确定的主题
LDA模型的介绍
模型是一个词袋模型,其认为文章是由一组词构成的一个集合.不关心文章中词的先后顺序.每一篇文章会由多个主题组成(主题分布θ),而对于每一个主题,又包括了各单词的出现概率φ
此时就可以获取某一个单词在文章中的出现概率:
p(w|d) = p(w|t) * p(t|d)
其中 p(w|t) 可以由分布 φ求得,p(t|d)可以由分布θ求得
模型的获取过程
最初提出LDA时:模型的更新使用的是EM算法,即先固定一个分布,对另一个分布进行更新,更新完成后再对另一个分布进行更新.(这也是sklearn中实现的方法)
对文档中的任意一个词,计算其在任意的topic下的p(w|t),然后就可以句topic对应的p(t|d),获取每一个opic下的p(w|d),取最大的那个,作为当前的词的topic.
具体过程为(来源于wiki百科)
吉布斯采样（Gibbs sampling）是统计学中用于马尔科夫蒙特卡洛（MCMC）的一种算法 .
蒙特卡洛和拉斯维加斯方法:
这是两类随机算法的统称
蒙特卡洛算法的思想在于每一次的采样结束后都会得到一个当前的最优算法,因此,采样越多获得的结果越接近最优解
拉斯维加斯方法则是有一个明确的求解目标,每一次的采样都会增大获得最优解的可能,因此这种方法可能会得不到解.
是基于多变量的条件概率分布,来近似抽取联合分布的样本序列.
也就是对多变量的条件分布进行连续采样的算法.
吉布斯采样的三种用途:
近似联合分布
部分变量的边缘分布
计算某一变量的期望值
具体的含义通过吉布斯采样的实现过程来解释.
在介绍采样的过程之前,需要先了解马尔科夫链
马尔科夫链
就是对于一个序列(包含先后顺序的样本集合),从上一个状态到下一个状态随机转移的过程.转移的过程是不具有记忆性的,即下一个时刻的状态只会由当前时刻的状态来决定,与历史的状态没有关系.因此就需要根据状态转移矩阵来获取下一时刻的状态.
而这个状态转移矩阵其实就是一个条件分布
采样过程
如果对全部的变量进行采样时,得到的就是所有变量的联合分布.
如果只有一部分变量需要变化,得到的是边缘分布
还可以对采样获得的所有样本求取某一变量的平均值来估计该变量的期望
首先基于先验概率 φ和θ,分别是 主题-词分布以及 文章-主题分布
联合分布概率是p(w,z),其中w是已知的,就是直接观测到的文本数据,而z则是隐含的主题分布.其实就是想要知道P(z|w)
对于每一个词的在其他所有词的词和话题的条件下,话题为k的概率为:
-i 表示去除词i之后剩余的其他所有词
文章中对于上述概率进行了详细的推导
最终得到:
其中n(k)m,¬i表示去除样本i后标记为话题k的词的数量
n(k)m,¬i则表示去除样本i后,词t的个数
首先随机初始化所有的词的topic(或者根据预设的分布来确定初始的topic)
重新扫描语料库,对每一个词重新试验吉布斯采样的格式获取他的topic
重复采样过程直到收敛
最终统计样本的topic-word矩阵以及 article-topic矩阵
Parameter estimation for text analysis
import re
import time
import jieba
import jieba.posseg as pseg
import numpy as np
import numpy
from sklearn.externals import joblib
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.decomposition import LatentDirichletAllocation
from sklearn.model_selection import GridSearchCV
def read_data(fname):
'''
输入的数据是按照:
title content answer tag tag1
存储的,在训练lda时,只需要使用title和answer.
:return: list 每一项是一个问题以及其所有的答案,使用.分隔
'''
dic = {}
with open(fname,'r',encoding='utf8') as f:
for i in f:
lst = i.strip().split('\t')
if lst[0] not in dic:
dic[lst[0]] = re.sub('[。？！，、….?!:]$','',lst[2])
else:
dic[lst[0]] += '.' + re.sub('[。？！，、….?!:]$','',lst[2])
corpus = [re.sub('[。？！，、….?!:]$','',key) + '.' + dic[key] for key in dic]
return corpus
def jieba_cut(corpus,cut_file,stopwords_file):
#结巴分词,并去除停用词(停用词为网上找到的中文停用词库),最后存储在cut_file中
stopwords = []
with open(stopwords_file,'r',encoding='utf8') as f:
for i in f:
stopwords.append(i.strip())
corpus_cut = []
n = 0
for s in corpus:
s_cut = [w for w in jieba.cut(s) if w not in stopwords]
corpus_cut.append(' '.join(s_cut))
n += 1
if n % 10000 == 0:
print(n)
f1 = open(cut_file, 'a', encoding='utf8')
for i in corpus_cut:
f1.write(i + '\n')
f1.close()
return corpus_cut
def vec_model(cut_file):
with open(cut_file, 'r', encoding='utf8') as f:
corpus_cut = [i.strip() for i in f.readlines()]
tf_vectorizer = CountVectorizer(max_df=0.95,min_df=2,stop_words='english')
x = tf_vectorizer.fit_transform(corpus_cut)
joblib.dump(tf_vectorizer,tf_ModelPath )
return x,tf_vectorizer,corpus_cut
def read_vec_model(cut_file,tf_ModelPath):
# 直接加载模型
with open(cut_file, 'r', encoding='utf8') as f:
corpus_cut = [i.strip() for i in f.readlines()]
tf_vectorizer = joblib.load(tf_ModelPath)
x = tf_vectorizer.fit_transform(corpus_cut)
return x,tf_vectorizer,corpus_cut
def train(vec_data,tf_model,n_topics = 14,max_iter = 10,learning_method= 'batch'):
'''
训练lda模型并存储
:param vec_data:
:param n_topics:
:param max_iter:
:param learning_method:
:return: 返回最终的lda模型
'''
lda = LatentDirichletAllocation(n_topics=n_topics,max_iter=max_iter,learning_method=learning_method,max_doc_update_iter=5)
print('train')
a = time.time()
lda.fit(vec_data)
print(time.time() - a)
n_top_words = 20
tf_feature_names = tf_model.get_feature_names()
print_top_words(lda, tf_feature_names, n_top_words)
joblib.dump(lda, lda_ModelPath)
return lda
def print_top_words(model, feature_names, n_top_words):
#打印每个主题下权重较高的term
for topic_idx, topic in enumerate(model.components_):
print( "Topic #%d:" % topic_idx)
print( " ".join([feature_names[i]for i in topic.argsort()[:-n_top_words - 1:-1]]))
print()
print(model.components_)
def grid_search(vec_data,tf_vectorizer,parameters):
GridSearchCV
def jieba_cut_transform(file):
#读取结巴分词的数据,然后将其中的英文和字母全部去除
with open(file,'r',encoding='utf8') as f:
ret_list = f.readlines()
f1 = open('./jieba_cut_all_drop.txt','a',encoding='utf8')
lst = []
for i in [[re.sub(r'[^\u4e00-\u9fa5]+','',j) for j in i.strip().split(' ') if re.sub(r'[a-zA-Z0-9]+','',j)] for i in ret_list]:
f1.write(' '.join(i) + '\n')
f1.close()
调用函数训练模型
lda_ModelPath = './lda_model2_all'
tf_ModelPath = './tf_model1_all'
fname = '../train_data'
cut_file = './jieba_cut_all.txt'
stopwords_file = './stopwords'
#读取数据
corpus = read_data(fname)
#分词
cut_data = jieba_cut(corpus, cut_file, stopwords_file)
#训练词向量模型
vec_data,tf_vectorizer,cut_data = vec_model(cut_file)
##这是直接读取数据的函数
#vec_data,tf_vectorizer,cut_data = read_vec_model(cut_file, tf_ModelPath)
lda = train(vec_data,tf_vectorizer,max_iter=50)
# 计算困惑度
# lda = joblib.load(lda_ModelPath)
p = lda.perplexity(vec_data)
print(p)
test_data = vec_data[:10]
ret = lda.transform(vec_data[:10])
print(ret)
print(ret.argmax(1))
print(cut_data[:10])
print(len(cut_data[:10]))
功能说明:
_check_params:
检查传入的参数是否满足要求
_init_latent_vars:
将输入的参数转换为训练时需要使用的参数
_e_step
em算法中的e步
_em_step
使用em算法迭代
_check_non_neg_array:
检查是否包含负数
partial_fit
Online VB with Mini-Batch update
fit
训练模型
transform
使用模型对数据进行预测
_approx_bound
估算log-likelihood
score
使用log-likelihood作为评价指标的分数
perplexity
计算混淆度指标
class LatentDirichletAllocation(BaseEstimator, TransformerMixin):
"""Latent Dirichlet Allocation with online variational Bayes algorithm
.. versionadded:: 0.17
Read more in the :ref:`User Guide <LatentDirichletAllocation>`.
Parameters
----------
n_topics : int, optional (default=10)
话题的数量。
doc_topic_prior : float, optional (default=None)
文档主题分布的先验‘theta’。如果值为None，
默认为“1 / n_topics”。
在文献中，这被称为“阿尔法”。
topic_word_prior : float, optional (default=None)
主题词分布的先验' beta '。如果该值为None，则为默认值
' 1 / n_topics '。
在文献中，这被称为“eta”。
learning_method : 'batch' | 'online', default='online'
方法用于更新' _component '。只用于' fit '方法。
一般来说，如果数据量很大，在线更新就会比批量更新快。
默认的学习方法将在0.20版本中更改为“batch”。
Valid options::
'batch': 分批变分贝叶斯法。使用所有的培训数据
每个EM更新。
旧的‘components_’将在每次迭代中被覆盖。
'online': 在线变分贝叶斯法。在每个EM更新中使用
更新“components_”的小批培训数据
变量增量。学习速率由
' learning_decay ' '和' learning_offset ' '参数。
learning_decay : float, optional (default=0.7)
它是在线学习方法中控制学习速率的一个参数
。值应该设置在(0.5,1.0)之间以保证
渐近收敛。当值为0.0且batch_size为n_samples时
，更新方法与批量学习相同。在文学，这叫做kappa。
learning_offset : float, optional (default=10.)
一个(正的)参数，它降低了在线学习中早期迭代的权重。它应该大于1.0。在文献中，这称为tau_0。
max_iter : integer, optional (default=10)
最大迭代次数。
total_samples : int, optional (default=1e6)
文件总数。只在' partial_fit '方法中使用。
batch_size : int, optional (default=128)
在每个EM迭代中使用的文档数量。仅用于在线学习。
evaluate_every : int optional (default=0)
在每个EM迭代中使用的文档数量。仅用于在线学习。
只用于' fit '方法。将其设为0或负数，则根本不能评价训练的复杂程度。对perplexity的评估可以帮助你在训练过程中检查收敛性，但是也会增加训练的总时间。在每次迭代中评估perplexity可能会增加两倍的训练时间。
perp_tol : float, optional (default=1e-1)
批量学习中的困惑容忍。仅当' ' evaluate_every ' '大于0时使用。
mean_change_tol : float, optional (default=1e-3)
停止E-step中更新文档主题分布的容忍。
max_doc_update_iter : int (default=100)
E-step中更新文档主题分布的最大迭代次数。
n_jobs : int, optional (default=1)
要在E-step中使用的作业数。如果是-1，则使用所有cpu。对于-1以下的' ' n_jobs ' '，则使用(n_cpu + 1 + n_jobs)。
verbose : int, optional (default=0)
冗长的水平。
random_state : int or RandomState instance or None, optional (default=None)
伪随机数发生器种子控制。
Attributes
----------
components_ : array, [n_topics, n_features]
主题字分布。components_[i, j] '代表主题i中的单词j。
n_batch_iter_ : int
EM步骤的迭代次数。
n_iter_ : int
通过数据集的次数。
"""
def partial_fit(self, X, y=None):
"""Online VB with Mini-Batch update.
Parameters
----------
X : array-like or sparse matrix, shape=(n_samples, n_features)
Document word matrix.
Returns
-------
self
"""
self._check_params()
X = self._check_non_neg_array(X,
"LatentDirichletAllocation.partial_fit")
n_samples, n_features = X.shape
batch_size = self.batch_size
# initialize parameters or check
if not hasattr(self, 'components_'):
self._init_latent_vars(n_features)
if n_features != self.components_.shape[1]:
raise ValueError(
"The provided data has %d dimensions while "
"the model was trained with feature size %d." %
(n_features, self.components_.shape[1]))
n_jobs = _get_n_jobs(self.n_jobs)
with Parallel(n_jobs=n_jobs, verbose=max(0, self.verbose - 1)) as parallel:
for idx_slice in gen_batches(n_samples, batch_size):
self._em_step(X[idx_slice, :],
total_samples=self.total_samples,
batch_update=False,
parallel=parallel)
return self
def fit(self, X, y=None):
"""Learn model for the data X with variational Bayes method.
When `learning_method` is 'online', use mini-batch update.
Otherwise, use batch update.
Parameters
----------
X : array-like or sparse matrix, shape=(n_samples, n_features)
Document word matrix.
Returns
-------
self
"""
#数据的准备
self._check_params()
X = self._check_non_neg_array(X, "LatentDirichletAllocation.fit")
n_samples, n_features = X.shape
max_iter = self.max_iter
evaluate_every = self.evaluate_every
learning_method = self.learning_method
if learning_method == None:
warnings.warn("The default value for 'learning_method' will be "
"changed from 'online' to 'batch' in the release 0.20. "
"This warning was introduced in 0.18.",
DeprecationWarning)
learning_method = 'online'
batch_size = self.batch_size
# initialize parameters
self._init_latent_vars(n_features)
# change to perplexity later
#上一次训练的困惑度
last_bound = None
n_jobs = _get_n_jobs(self.n_jobs)
with Parallel(n_jobs=n_jobs, verbose=max(0, self.verbose - 1)) as parallel:
for i in xrange(max_iter):
#在online模式下,划分为batch_size大小的数据集训练,batch_update为False
#在batch模式下batch_update为True
if learning_method == 'online':
for idx_slice in gen_batches(n_samples, batch_size):
self._em_step(X[idx_slice, :], total_samples=n_samples,
batch_update=False, parallel=parallel)
else:
# batch update
self._em_step(X, total_samples=n_samples,
batch_update=True, parallel=parallel)
# check perplexity,根据evaluate_every的设定值,来判断是否在训练的时候检查困惑度以及检查的频率,检查困惑度会花费大量的时间.
if evaluate_every > 0 and (i + 1) % evaluate_every == 0:
doc_topics_distr, _ = self._e_step(X, cal_sstats=False,
random_init=False,
parallel=parallel)
bound = self.perplexity(X, doc_topics_distr,
sub_sampling=False)
if self.verbose:
print('iteration: %d, perplexity: %.4f'
% (i + 1, bound))
if last_bound and abs(last_bound - bound) < self.perp_tol:
break
last_bound = bound
self.n_iter_ += 1
return self
使用到的函数:
def _check_non_neg_array(self, X, whom):
"""check X format
check X format and make sure no negative value in X.
Parameters
----------
X : array-like or sparse matrix
"""
X = check_array(X, accept_sparse='csr')
#检查是否包含负数
check_non_negative(X, whom)
"""
check_non_negative:
X = X.data if sp.issparse(X) else X
if (X < 0).any():
raise ValueError("Negative values in data passed to %s" % whom)
"""
return X
"""
Parallel:
并行处理数据的类,使用到的两个参数的作用为:
n_jobs:
并发运行作业的最大数量，例如后端=“多处理”时Python工作进程的数量或后端=“线程化”时线程池的大小。如果-1，则使用所有cpu。如果给定1，则根本不使用并行计算代码，这对于调试非常有用。对于小于-1的n_jobs，则使用(n_cpu + 1 + n_jobs)。因此，对于n_jobs = -2，将使用除一个cpu之外的所有cpu。
verbose:
冗余级别:如果非零，则打印进度消息。在50以上，输出被发送到stdout。消息的频率随冗长级别的增加而增加。如果大于10，则报告所有迭代。
"""
def gen_batches(n, batch_size):
"""Generator to create slices containing batch_size elements, from 0 to n.
The last slice may contain less than batch_size elements, when batch_size
does not divide n.
Examples
--------
>>> from sklearn.utils import gen_batches
>>> list(gen_batches(7, 3))
[slice(0, 3, None), slice(3, 6, None), slice(6, 7, None)]
>>> list(gen_batches(6, 3))
[slice(0, 3, None), slice(3, 6, None)]
>>> list(gen_batches(2, 3))
[slice(0, 2, None)]
"""
start = 0
for _ in range(int(n // batch_size)):
end = start + batch_size
yield slice(start, end)
start = end
if start < n:
yield slice(start, n)
def _em_step(self, X, total_samples, batch_update, parallel=None):
"""EM update for 1 iteration.
update `_component` by batch VB or online VB.
Parameters
----------
X : array-like or sparse matrix, shape=(n_samples, n_features)
Document word matrix.
total_samples : integer
Total umber of documents. It is only used when batch_update is `False`.
batch_update : boolean
Parameter that controls updating method.
`True` for batch learning, `False` for online learning.
parallel : joblib.Parallel
预先初始化的joblib.Parallel实例。
Returns
-------
doc_topic_distr : array, shape=(n_samples, n_topics)
非规范化文档主题分布.
"""
# E-step
#返回值是文章-话题分布:
_, suff_stats = self._e_step(X, cal_sstats=True, random_init=True,
parallel=parallel)
# M-step
#如果是batch模式,则
if batch_update:
#使用刚刚得到的文章话题分布更新话题词频的分布
self.components_ = self.topic_word_prior_ + suff_stats
else:
# online update
# In the literature, the weight is `rho`
#np.power第一个参数对第二个参数取幂
weight = np.power(self.learning_offset + self.n_batch_iter_,
-self.learning_decay)
doc_ratio = float(total_samples) / X.shape[0]
#上一次的结果变为1-w倍,当前的累加值系数为w
self.components_ *= (1 - weight)
self.components_ += (weight * (self.topic_word_prior_
+ doc_ratio * suff_stats))
# update `component_` related variables
self.exp_dirichlet_component_ = np.exp(
_dirichlet_expectation_2d(self.components_))
self.n_batch_iter_ += 1
return
E步:
def _e_step(self, X, cal_sstats, random_init, parallel=None):
"""E-step in EM update.
Parameters
----------
X : array-like or sparse matrix, shape=(n_samples, n_features)
Document word matrix.
cal_sstats : boolean
参数，指示是否计算足够的统计信息。当需要运行M-step时，将' ' cal_sstats ' '设置为True。
random_init : boolean
参数，指示是否在e -步骤中随机初始化文档主题分布。在训练步骤中将它设置为True。
parallel : joblib.Parallel (optional)
预先初始化的joblib.Parallel实例。
Returns
-------
(doc_topic_distr, suff_stats) :
`doc_topic_distr` 为每个文档的非规范化主题分布。在文献中，这叫做 `gamma`.
`suff_stats` 预计M-step有足够的统计数据。当' cal_sstats == False '时，它将为None。
"""
# Run e-step in parallel
random_state = self.random_state_ if random_init else None
# TODO: make Parallel._effective_n_jobs public instead?
# 获取计算所需的作业数。
n_jobs = _get_n_jobs(self.n_jobs)
if parallel is None:
#如果为空则创建一个parallel对象
parallel = Parallel(n_jobs=n_jobs, verbose=max(0, self.verbose - 1))
#使用parallel对象并行执行更新过程
#delayed是parallel中用于捕获函数参数的装饰器。
#_update_doc_distribution则是最终需要调用来更新的函数
#parallel的call方法需要的参数是要执行的函数的可迭代序列
results = parallel(
delayed(_update_doc_distribution)(X[idx_slice, :],
self.exp_dirichlet_component_,
self.doc_topic_prior_,
self.max_doc_update_iter,
self.mean_change_tol, cal_sstats,
random_state)
for idx_slice in gen_even_slices(X.shape[0], n_jobs))
# merge result
doc_topics, sstats_list = zip(*results)
doc_topic_distr = np.vstack(doc_topics)
if cal_sstats:
# This step finishes computing the sufficient statistics for the
# M-step.
suff_stats = np.zeros(self.components_.shape)
for sstats in sstats_list:
suff_stats += sstats
suff_stats *= self.exp_dirichlet_component_
else:
suff_stats = None
return (doc_topic_distr, suff_stats)
def _update_doc_distribution(X, exp_topic_word_distr, doc_topic_prior,
max_iters,
mean_change_tol, cal_sstats, random_state):
"""E-step: 更新文档主题分布.
Parameters
----------
X : array-like or sparse matrix, shape=(n_samples, n_features)
文档-词矩阵.
exp_topic_word_distr : dense matrix, shape=(n_topics, n_features)
对数主题词分布期望的指数值。在文献中，这是 `exp(E[log(beta)])`.
doc_topic_prior : float
文档主题分发的优先级 `theta`.
max_iters : int
Max number of iterations for updating document topic distribution in
the E-step.
mean_change_tol : float
停止E-setp中更新文档主题分布的容忍。
cal_sstats : boolean
指示是否计算足够统计量的参数。
Set `cal_sstats` to `True` when we need to run M-step.
random_state : RandomState instance or None
参数，指示如何初始化文档主题分布。
Set `random_state` to None will initialize document topic distribution
to a 常数.
Returns
-------
(doc_topic_distr, suff_stats) :
`doc_topic_distr` is unnormalized topic distribution for each document.
In the literature, this is `gamma`. we can calculate `E[log(theta)]`
from it.
`suff_stats` is expected sufficient statistics for the M-step.
When `cal_sstats == False`, this will be None.
"""
is_sparse_x = sp.issparse(X)
n_samples, n_features = X.shape
n_topics = exp_topic_word_distr.shape[0]
if random_state:
doc_topic_distr = random_state.gamma(100., 0.01, (n_samples, n_topics))
else:
doc_topic_distr = np.ones((n_samples, n_topics))
# In the literature, this is `exp(E[log(theta)])`
exp_doc_topic = np.exp(_dirichlet_expectation_2d(doc_topic_distr))
# diff on `component_` (only calculate it when `cal_diff` is True)
suff_stats = np.zeros(exp_topic_word_distr.shape) if cal_sstats else None
if is_sparse_x:
#是否是稀疏的
X_data = X.data
X_indices = X.indices
X_indptr = X.indptr
for idx_d in xrange(n_samples):
if is_sparse_x:
ids = X_indices[X_indptr[idx_d]:X_indptr[idx_d + 1]]
cnts = X_data[X_indptr[idx_d]:X_indptr[idx_d + 1]]
else:
ids = np.nonzero(X[idx_d, :])[0]
cnts = X[idx_d, ids]
doc_topic_d = doc_topic_distr[idx_d, :]
# The next one is a copy, since the inner loop overwrites it.
exp_doc_topic_d = exp_doc_topic[idx_d, :].copy()
exp_topic_word_d = exp_topic_word_distr[:, ids]
# Iterate between `doc_topic_d` and `norm_phi` until convergence
for _ in xrange(0, max_iters):
#迭代max_iters次
last_d = doc_topic_d
#迭代公式如下:
# The optimal phi_{dwk} is proportional to
# exp(E[log(theta_{dk})]) * exp(E[log(beta_{dw})]).
norm_phi = np.dot(exp_doc_topic_d, exp_topic_word_d) + EPS
doc_topic_d = (exp_doc_topic_d *
np.dot(cnts / norm_phi, exp_topic_word_d.T))
# Note: adds doc_topic_prior to doc_topic_d, in-place.
_dirichlet_expectation_1d(doc_topic_d, doc_topic_prior,
exp_doc_topic_d)
if mean_change(last_d, doc_topic_d) < mean_change_tol:
#如果两次的变化值小于设定的阈值,停止迭代
break
doc_topic_distr[idx_d, :] = doc_topic_d
# Contribution of document d to the expected sufficient
# statistics for the M step.
if cal_sstats:
norm_phi = np.dot(exp_doc_topic_d, exp_topic_word_d) + EPS
#np.outer两个向量的外积
suff_stats[:, ids] += np.outer(exp_doc_topic_d, cnts / norm_phi)
return (doc_topic_distr, suff_stats)
def transform(self, X):
"""Transform data X according to the fitted model.
Parameters
----------
X : array-like or sparse matrix, shape=(n_samples, n_features)
Document word matrix.
Returns
-------
doc_topic_distr : shape=(n_samples, n_topics)
Document topic distribution for X.
"""
if not hasattr(self, 'components_'):
#判断是否存在两个分布矩阵
raise NotFittedError("no 'components_' attribute in model."
" Please fit model first.")
# make sure feature size is the same in fitted model and in X
#对传入数据的特征数量
X = self._check_non_neg_array(X, "LatentDirichletAllocation.transform")
n_samples, n_features = X.shape
if n_features != self.components_.shape[1]:
raise ValueError(
"The provided data has %d dimensions while "
"the model was trained with feature size %d." %
(n_features, self.components_.shape[1]))
doc_topic_distr, _ = self._e_step(X, cal_sstats=False,random_init=False)
# normalize doc_topic_distr
doc_topic_distr /= doc_topic_distr.sum(axis=1)[:, np.newaxis]
return doc_topic_distr
def _approx_bound(self, X, doc_topic_distr, sub_sampling):
"""Estimate the variational bound.
仅使用作为x传递的文档估计“所有文档”的变分界，因为每个单词的log-likelihood不能直接计算，所以我们使用这个界来估计。
Parameters
----------
X : array-like or sparse matrix, shape=(n_samples, n_features)
Document word matrix.
doc_topic_distr : array, shape=(n_samples, n_topics)
Document topic distribution. In the literature, this is called
gamma.
sub_sampling : boolean, optional, (default=False)
Compensate for subsampling of documents.
It is used in calculate bound in online learning.
Returns
-------
score : float
"""
def _loglikelihood(prior, distr, dirichlet_distr, size):
# calculate log-likelihood
score = np.sum((prior - distr) * dirichlet_distr)
score += np.sum(gammaln(distr) - gammaln(prior))
score += np.sum(gammaln(prior * size) - gammaln(np.sum(distr, 1)))
return score
is_sparse_x = sp.issparse(X)
n_samples, n_topics = doc_topic_distr.shape
n_features = self.components_.shape[1]
score = 0
dirichlet_doc_topic = _dirichlet_expectation_2d(doc_topic_distr)
dirichlet_component_ = _dirichlet_expectation_2d(self.components_)
doc_topic_prior = self.doc_topic_prior_
topic_word_prior = self.topic_word_prior_
if is_sparse_x:
X_data = X.data
X_indices = X.indices
X_indptr = X.indptr
# E[log p(docs | theta, beta)]
for idx_d in xrange(0, n_samples):
if is_sparse_x:
ids = X_indices[X_indptr[idx_d]:X_indptr[idx_d + 1]]
cnts = X_data[X_indptr[idx_d]:X_indptr[idx_d + 1]]
else:
ids = np.nonzero(X[idx_d, :])[0]
cnts = X[idx_d, ids]
temp = (dirichlet_doc_topic[idx_d, :, np.newaxis]
+ dirichlet_component_[:, ids])
norm_phi = logsumexp(temp)
score += np.dot(cnts, norm_phi)
# compute E[log p(theta | alpha) - log q(theta | gamma)]
score += _loglikelihood(doc_topic_prior, doc_topic_distr,
dirichlet_doc_topic, self.n_topics)
# Compensate for the subsampling of the population of documents
if sub_sampling:
doc_ratio = float(self.total_samples) / n_samples
score *= doc_ratio
# E[log p(beta | eta) - log q (beta | lambda)]
score += _loglikelihood(topic_word_prior, self.components_,
dirichlet_component_, n_features)
return score
################################################################################
def score(self, X, y=None):
"""Calculate approximate log-likelihood as score.
Parameters
----------
X : array-like or sparse matrix, shape=(n_samples, n_features)
Document word matrix.
Returns
-------
score : float
Use approximate bound as score.
"""
X = self._check_non_neg_array(X, "LatentDirichletAllocation.score")
doc_topic_distr = self.transform(X)
score = self._approx_bound(X, doc_topic_distr, sub_sampling=False)
return score
################################################################################
def perplexity(self, X, doc_topic_distr=None, sub_sampling=False):
"""Calculate approximate perplexity for data X.
Perplexity is defined as exp(-1. * log-likelihood per word)
Parameters
----------
X : array-like or sparse matrix, [n_samples, n_features]
Document word matrix.
doc_topic_distr : None or array, shape=(n_samples, n_topics)
Document topic distribution.
If it is None, it will be generated by applying transform on X.
Returns
-------
score : float
Perplexity score.
"""
if not hasattr(self, 'components_'):
raise NotFittedError("no 'components_' attribute in model."
" Please fit model first.")
X = self._check_non_neg_array(X,
"LatentDirichletAllocation.perplexity")
if doc_topic_distr is None:
doc_topic_distr = self.transform(X)
else:
n_samples, n_topics = doc_topic_distr.shape
if n_samples != X.shape[0]:
raise ValueError("Number of samples in X and doc_topic_distr"
" do not match.")
if n_topics != self.n_topics:
raise ValueError("Number of topics does not match.")
current_samples = X.shape[0]
bound = self._approx_bound(X, doc_topic_distr, sub_sampling)
if sub_sampling:
word_cnt = X.sum() * (float(self.total_samples) / current_samples)
else:
word_cnt = X.sum()
perword_bound = bound / word_cnt
return np.exp(-1.0 * perword_bound)
根据源码可以看到其包含的参数如下:
n_topics : 话题的数量。
doc_topic_prior : 文档主题分布的先验
topic_word_prior : 主题词分布的先验
learning_method : 'batch' | 'online', default='online'
'batch': 分批变分贝叶斯法。使用所有的培训数据更新每个EM。
旧的‘components_’将在每次迭代中被覆盖。
'online': 在线变分贝叶斯法。在每个EM中使用小批量数据更新“components_”的变量增量。学习速率由 learning_decay和learning_offset。
learning_decay :(default=0.7)在线学习方法中控制学习速率的一个参数。值应该设置在(0.5,1.0)之间以保证渐近收敛。当值为0.0且batch_size为n_samples时，更新方法与批量学习相同。
learning_offset : (default=10.)一个(正的)参数，它降低了在线学习中早期迭代的权重。它应该大于1.0。
max_iter : (default=10)最大迭代次数。
total_samples :(default=1e6)文件总数。只在' partial_fit '方法中使用。
batch_size : (default=128)在每个EM迭代中使用的文档数量。仅用于在线学习。
evaluate_every : (default=0)在每个EM迭代中使用的文档数量。仅用于在线学习。
只用于' fit '方法。将其设为0或负数，则根本不能评价训练的复杂程度。对perplexity的评估可以帮助你在训练过程中检查收敛性，但是也会增加训练的总时间。在每次迭代中评估perplexity可能会增加两倍的训练时间。(设置这个参数后,会在循环的计数器i对evaluate_every取余为0时,输出一次混淆度)
perp_tol : float, optional (default=1e-1)
批量学习中的困惑容忍。仅当evaluate_every大于0时使用。
mean_change_tol :(default=1e-3)
停止E-step中更新文档主题分布的容忍。
max_doc_update_iter : (default=100)
E-step中更新文档主题分布的最大迭代次数。
n_jobs : int, optional (default=1)
要在E-step中使用的作业数。如果是-1，则使用所有cpu。对于-1以下的' ' n_jobs ' '，则使用(n_cpu + 1 + n_jobs)。
verbose : int, optional (default=0)
冗长的水平。
random_state : (default=None) 伪随机数发生器种子控制。
Attributes
----------
components_ : array, [n_topics, n_features]
主题字分布。components_[i, j] 代表主题i中的单词j。
n_batch_iter_ : int
EM步骤的迭代次数。
n_iter_ : int
通过数据集的次数。
因此在训练时主要需要调节的参数就是:
1.n_topics这就相当于聚类的个数
2.learning_method 这个是用来控制训练的模式的,如果选择了online模式,还需要对其他的一些参数进行调节,来控制训练的学习率等参数
3.最大迭代次数,包括每一次em算法中的停止迭代的条件和迭代次数以及对于数据集的迭代次数
4.n_jobs是调节并行迭代的参数
5.如果想在训练的过程中观察模型的指标的变化情况,需要对evaluate_every进行设置,但是不会影响模型的结果,训练的时间反而会大幅度的增加.