Create a gist now

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Semi-supervised Naive Bayes
# -*- coding: utf-8 -*-
# Copyright (C) 2010 Mathieu Blondel
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Implementation details are described at
# http://www.mblondel.org/journal/2010/06/21/semi-supervised-naive-bayes-in-python/
import numpy as np
"""
Implementation of Naive Bayes trained by EM for semi-supervised text
classification, as described in
"Semi-Supervised Text Classification Using EM", by Nigam et al.
Notation:
w: word
d: document
c: class
V: vocabulary size
X: number of documents
M: number of classes
"""
def softmax(loga, k=-np.inf, out=None):
"""
Compute the sotfmax function (normalized exponentials) without underflow.
exp^a_i / \sum_j exp^a_j
"""
if out is None: out = np.empty_like(loga).astype(np.double)
m = np.max(loga)
logam = loga - m
sup = logam > k
inf = np.logical_not(sup)
out[sup] = np.exp(logam[sup])
out[inf] = 0.0
out /= np.sum(out)
return out
def logsum(loga, k=-np.inf):
"""
Compute a sum of logs without underflow.
\log \sum_c e^{b_c} = log [(\sum_c e^{b_c}) e^{-B}e^B]
= log [(\sum_c e^{b_c-B}) e^B]
= [log(\sum_c e^{b_c-B}) + B
where B = max_c b_c
"""
B = np.max(loga)
logaB = aB = loga - B
sup = logaB > k
inf = np.logical_not(sup)
aB[sup] = np.exp(logaB[sup])
aB[inf] = 0.0
return (np.log(np.sum(aB)) + B)
def loglikelihood(td, delta, tdu, p_w_c_log, p_c_log):
V, Xl = td.shape
V_, Xu = tdu.shape
Xl_, M = delta.shape
lik = 0.0
## labeled
# log P(x|c)
p_x_c_log = np.zeros((Xl,M), np.double)
for w,d in zip(*td.nonzero()):
p_x_c_log[d,:] += p_w_c_log[w,:] * td[w,d]
# add log P(c) + lop P(x|c) if x has label c
for d,c in zip(*delta.nonzero()):
lik += p_c_log[c] + p_x_c_log[d,c]
## unlabelled
# log P(x|c)
p_x_c_log = np.zeros((Xu,M), np.double)
for w,d in zip(*tdu.nonzero()):
p_x_c_log[d,:] += p_w_c_log[w,:] * tdu[w,d]
# add log P(c)
p_x_c_log += p_c_log[np.newaxis,:]
for d in range(Xu):
lik += logsum(p_x_c_log[d,:], k=-10)
return lik
def normalize_p_c(p_c):
M = len(p_c)
denom = M + np.sum(p_c)
p_c += 1.0
p_c /= denom
def normalize_p_w_c(p_w_c):
V, X = p_w_c.shape
denoms = V + np.sum(p_w_c, axis=0)
p_w_c += 1.0
p_w_c /= denoms[np.newaxis,:]
class SemiNB(object):
def __init__(self, model=None):
"""
model: a model, as returned by get_model() or train().
"""
self.p_w_c = None
self.p_c = None
if model is not None: self.set_model(model)
self.debug = False
def train(self, td, delta, normalize=True, sparse=True):
"""
td: term-document matrix V x X
delta: X x M matrix
where delta(d,c) = 1 if document d belongs to class c
"""
X_, M = delta.shape
V, X = td.shape
assert(X_ == X)
# P(c)
self.p_c = np.sum(delta, axis=0)
# P(w|c)
self.p_w_c = np.zeros((V,M), dtype=np.double)
if sparse:
# faster when delta is sparse
# select indices of documents that have class c
for d,c in zip(*delta.nonzero()):
# select indices of terms that are non-zero
for w in np.flatnonzero(td[:,d]):
self.p_w_c[w,c] += td[w,d] * delta[d,c]
else:
# faster when delta is non-sparse
for w,d in zip(*td.nonzero()):
self.p_w_c[w,:] += td[w,d] * delta[d,:]
if normalize:
normalize_p_c(self.p_c)
normalize_p_w_c(self.p_w_c)
return self.get_model()
def train_semi(self, td, delta, tdu, maxiter=50, eps=0.01):
"""
td: V x X term document matrix
delta: X x M label matrix
tdu: V x Xu term document matrix (unlabeled)
maxiter: maximum number of iterations
eps: stop if no more progress than esp
"""
X_, M = delta.shape
V, X = td.shape
assert(X_ == X)
# compute counts for labeled data once for all
self.train(td, delta, normalize=False)
p_c_l = np.array(self.p_c, copy=True)
p_w_c_l = np.array(self.p_w_c, copy=True)
# normalize to get initial classifier
normalize_p_c(self.p_c)
normalize_p_w_c(self.p_w_c)
lik = loglikelihood(td, delta, tdu, np.log(self.p_w_c), np.log(self.p_c))
for iteration in range(1, maxiter+1):
# E-step: find the probabilistic labels for unlabeled data
delta_u = self.predict_proba_all(tdu)
# M-step: train classifier with the union of
# labeled and unlabeled data
self.train(tdu, delta_u, normalize=False, sparse=False)
self.p_c += p_c_l
self.p_w_c += p_w_c_l
normalize_p_c(self.p_c)
normalize_p_w_c(self.p_w_c)
lik_new = loglikelihood(td, delta, tdu,
np.log(self.p_w_c), np.log(self.p_c))
lik_diff = lik_new - lik
assert(lik_diff >= -1e-10)
lik = lik_new
if lik_diff < eps:
print "No more progress, stopping EM at iteration", iteration
break
if self.debug:
print "Iteration", iteration
print "L += %f" % lik_diff
return self.get_model()
def p_x_c_log_all(self, td):
M = len(self.p_c)
V, X = td.shape
p_x_c_log = np.zeros((X,M), np.double)
p_w_c_log = np.log(self.p_w_c)
# log P(x|c)
for w,d in zip(*td.nonzero()):
p_x_c_log[d,:] += p_w_c_log[w,:] * td[w,d]
return p_x_c_log
def predict_proba(self, x):
"""
x: a V array representing a document
Compute a M array containing P(c|x).
"""
return self.predict_proba_all(x[:,np.newaxis])[0,:]
def predict_proba_all(self, td):
"""
td: V x X term document matrix
Compute an X x M matrix of P(c|x) for all x in td.
"""
V, X = td.shape
# log P(x|c)
p_x_c_log = self.p_x_c_log_all(td)
# add log P(c)
p_x_c_log += np.log(self.p_c)[np.newaxis,:]
# sotfmax(log P(x|c) + log P(c)) = P(c|x)
for d in range(X):
softmax(p_x_c_log[d,:], k=-10, out=p_x_c_log[d,:])
return p_x_c_log
def predict(self, x):
"""
x: a V array representing a document
Compute the predicted class index.
"""
return self.predict_all(x[:,np.newaxis])[0]
def predict_all(self, td):
"""
td: V x X term document matrix
Compute a X array containing predicted class indices.
Note: the main difference with predict_proba_all is that the
normalization is not necessary, as we are only interested in the most
likely class.
"""
# log P(x|c)
p_x_c_log = self.p_x_c_log_all(td)
# add log P(c)
p_x_c_log += np.log(self.p_c)[np.newaxis,:]
return p_x_c_log.argmax(axis=1)
def get_model(self):
return (self.p_w_c, self.p_c)
def set_model(self, model):
self.p_w_c, self.p_c = model
@jinyyy666

This comment has been minimized.

Show comment
Hide comment
@jinyyy666

jinyyy666 Apr 4, 2017

Nice code!
In Line:113, should it be "V,M = p_w_c.shape" instead of "V, X = p_w_c.shape" for the sake of clarity although the "X" is not used here? : )

jinyyy666 commented Apr 4, 2017

Nice code!
In Line:113, should it be "V,M = p_w_c.shape" instead of "V, X = p_w_c.shape" for the sake of clarity although the "X" is not used here? : )

@amitsasane007

This comment has been minimized.

Show comment
Hide comment
@amitsasane007

amitsasane007 May 10, 2017

Please add description to code

amitsasane007 commented May 10, 2017

Please add description to code

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment