mvpa.clfs.gnb

33 """Gaussian Naive Bayes `Classifier`. 34 35 GNB is a probabilistic classifier relying on Bayes rule to 36 estimate posterior probabilities of labels given the data. Naive 37 assumption in it is an independence of the features, which allows 38 to combine per-feature likelihoods by a simple product across 39 likelihoods of"independent" features. 40 See http://en.wikipedia.org/wiki/Naive_bayes for more information. 41 42 Provided here implementation is "naive" on its own -- various 43 aspects could be improved, but has its own advantages: 44 45 - implementation is simple and straightforward 46 - no data copying while considering samples of specific class 47 - provides alternative ways to assess prior distribution of the 48 classes in the case of unbalanced sets of samples (see parameter 49 `prior`) 50 - makes use of NumPy broadcasting mechanism, so should be 51 relatively efficient 52 - should work for any dimensionality of samples 53 54 GNB is listed both as linear and non-linear classifier, since 55 specifics of separating boundary depends on the data and/or 56 parameters: linear separation is achieved whenever samples are 57 balanced (or prior='uniform') and features have the same variance 58 across different classes (i.e. if common_variance=True to enforce 59 this). 60 61 Whenever decisions are made based on log-probabilities (parameter 62 logprob=True, which is the default), then state variable `values` 63 if enabled would also contain log-probabilities. Also mention 64 that normalization by the evidence (P(data)) is disabled by 65 default since it has no impact per se on classification decision. 66 You might like set parameter normalize to True if you want to 67 access properly scaled probabilities in `values` state variable. 68 """ 69 # XXX decide when should we set corresponding internal, 70 # since it depends actually on the data -- no clear way, 71 # so set both linear and non-linear 72 _clf_internals = [ 'gnb', 'linear', 'non-linear', 73 'binary', 'multiclass' ] 74 75 common_variance = Parameter(False, allowedtype='bool', 76 doc="""Use the same variance across all classes.""") 77 prior = Parameter('laplacian_smoothing', 78 allowedtype='basestring', 79 choices=["laplacian_smoothing", "uniform", "ratio"], 80 doc="""How to compute prior distribution.""") 81 82 logprob = Parameter(True, allowedtype='bool', 83 doc="""Operate on log probabilities. Preferable to avoid unneeded 84 exponentiation and loose precision. 85 If set, logprobs are stored in `values`""") 86 normalize = Parameter(False, allowedtype='bool', 87 doc="""Normalize (log)prob by P(data). Requires probabilities thus 88 for `logprob` case would require exponentiation of 'logprob's, thus 89 disabled by default since does not impact classification output. 90 """) 91

92 - def __init__(self, **kwargs):

93 """Initialize an GNB classifier. 94 """ 95 96 # init base class first 97 Classifier.__init__(self, **kwargs) 98 99 # pylint friendly initializations 100 self.means = None 101 """Means of features per class""" 102 self.variances = None 103 """Variances per class, but "vars" is taken ;)""" 104 self.ulabels = None 105 """Labels classifier was trained on""" 106 self.priors = None 107 """Class probabilities""" 108 109 # Define internal state of classifier 110 self._norm_weight = None

111

112 - def _train(self, dataset):

113 """Train the classifier using `dataset` (`Dataset`). 114 """ 115 params = self.params 116 117 # get the dataset information into easy vars 118 X = dataset.samples 119 labels = dataset.labels 120 self.ulabels = ulabels = dataset.uniquelabels 121 nlabels = len(ulabels) 122 #params = self.params # for quicker access 123 label2index = dict((l, il) for il, l in enumerate(ulabels)) 124 125 # set the feature dimensions 126 nsamples = len(X) 127 s_shape = X.shape[1:] # shape of a single sample 128 129 self.means = means = \ 130 N.zeros((nlabels, ) + s_shape) 131 self.variances = variances = \ 132 N.zeros((nlabels, ) + s_shape) 133 # degenerate dimension are added for easy broadcasting later on 134 nsamples_per_class = N.zeros((nlabels,) + (1,)*len(s_shape)) 135 136 # Estimate means and number of samples per each label 137 for s, l in zip(X, labels): 138 il = label2index[l] # index of the label 139 nsamples_per_class[il] += 1 140 means[il] += s 141 142 # helped function - squash all dimensions but 1 143 squash = lambda x: N.atleast_1d(x.squeeze()) 144 ## Actually compute the means 145 non0labels = (squash(nsamples_per_class) != 0) 146 means[non0labels] /= nsamples_per_class[non0labels] 147 148 # Estimate variances 149 # better loop than repmat! ;) 150 for s, l in zip(X, labels): 151 il = label2index[l] # index of the label 152 variances[il] += (s - means[il])**2 153 154 ## Actually compute the variances 155 if params.common_variance: 156 # we need to get global std 157 cvar = N.sum(variances, axis=0)/nsamples # sum across labels 158 # broadcast the same variance across labels 159 variances[:] = cvar 160 else: 161 variances[non0labels] /= nsamples_per_class[non0labels] 162 163 # Store prior probabilities 164 prior = params.prior 165 if prior == 'uniform': 166 self.priors = N.ones((nlabels,))/nlabels 167 elif prior == 'laplacian_smoothing': 168 self.priors = (1+squash(nsamples_per_class)) \ 169 / (float(nsamples) + nlabels) 170 elif prior == 'ratio': 171 self.priors = squash(nsamples_per_class) / float(nsamples) 172 else: 173 raise ValueError( 174 "No idea on how to handle '%s' way to compute priors" 175 % params.prior) 176 177 # Precompute and store weighting coefficient for Gaussian 178 if params.logprob: 179 # it would be added to exponent 180 self._norm_weight = -0.5 * N.log(2*N.pi*variances) 181 else: 182 self._norm_weight = 1.0/N.sqrt(2*N.pi*variances) 183 184 if __debug__ and 'GNB' in debug.active: 185 debug('GNB', "training finished on data.shape=%s " % (X.shape, ) 186 + "min:max(data)=%f:%f" % (N.min(X), N.max(X)))

187 188

189 - def untrain(self):

190 """Untrain classifier and reset all learnt params 191 """ 192 self.means = None 193 self.variances = None 194 self.ulabels = None 195 self.priors = None 196 super(GNB, self).untrain()

197 198

199 - def _predict(self, data):

200 """Predict the output for the provided data. 201 """ 202 params = self.params 203 # argument of exponentiation 204 scaled_distances = \ 205 -0.5 * (((data - self.means[:, N.newaxis, ...])**2) \ 206 / self.variances[:, N.newaxis, ...]) 207 if params.logprob: 208 # if self.params.common_variance: 209 # XXX YOH: 210 # For decision there is no need to actually compute 211 # properly scaled p, ie 1/sqrt(2pi * sigma_i) could be 212 # simply discarded since it is common across features AND 213 # classes 214 # For completeness -- computing everything now even in logprob 215 lprob_csfs = self._norm_weight[:, N.newaxis, ...] + scaled_distances 216 217 # XXX for now just cut/paste with different operators, but 218 # could just bind them and reuse in the same equations 219 # Naive part -- just a product of probabilities across features 220 ## First we need to reshape to get class x samples x features 221 lprob_csf = lprob_csfs.reshape( 222 lprob_csfs.shape[:2] + (-1,)) 223 ## Now -- sum across features 224 lprob_cs = lprob_csf.sum(axis=2) 225 226 # Incorporate class probabilities: 227 prob_cs_cp = lprob_cs + N.log(self.priors[:, N.newaxis]) 228 229 else: 230 # Just a regular Normal distribution with per 231 # feature/class mean and variances 232 prob_csfs = \ 233 self._norm_weight[:, N.newaxis, ...] * N.exp(scaled_distances) 234 235 # Naive part -- just a product of probabilities across features 236 ## First we need to reshape to get class x samples x features 237 prob_csf = prob_csfs.reshape( 238 prob_csfs.shape[:2] + (-1,)) 239 ## Now -- product across features 240 prob_cs = prob_csf.prod(axis=2) 241 242 # Incorporate class probabilities: 243 prob_cs_cp = prob_cs * self.priors[:, N.newaxis] 244 245 # Normalize by evidence P(data) 246 if params.normalize: 247 if params.logprob: 248 prob_cs_cp_real = N.exp(prob_cs_cp) 249 else: 250 prob_cs_cp_real = prob_cs_cp 251 prob_s_cp_marginals = N.sum(prob_cs_cp_real, axis=0) 252 if params.logprob: 253 prob_cs_cp -= N.log(prob_s_cp_marginals) 254 else: 255 prob_cs_cp /= prob_s_cp_marginals 256 257 # Take the class with maximal (log)probability 258 winners = prob_cs_cp.argmax(axis=0) 259 predictions = [self.ulabels[c] for c in winners] 260 261 262 self.values = prob_cs_cp.T # set to the probabilities per class 263 264 if __debug__ and 'GNB' in debug.active: 265 debug('GNB', "predict on data.shape=%s min:max(data)=%f:%f " % 266 (data.shape, N.min(data), N.max(data))) 267 268 return predictions

Source Code for Module mvpa.clfs.gnb