Re: [mvpa-toolbox] Revisiting the issue of Gaussian Naive Bayes classifier- new and Final one.... I hope!

[Francisco Pereira]

I'm not entirely familiar with the way GNB is implemented in the toolbox, but what you describe sounds like what people generally do: they work in log-space, without using the normalization constant. Essentially, you want to compare (x is a test example, A|B are classes)

P(A|x) = P(x|A) P(A) / P(x)

with

P(B|x) = P(x|B) P(B) / P(x)

P(x) can be removed, as it is the same for both P(A|x) and P(B|x). As log is monotonic,
you can just compare

log P(x|A) P(A) = log P(x|A) + log P(A)

and

log P(x|B) P(B) = log P(x|B) + log P(B)

which seems to be what you are doing.

cheers,
Francisco

On Mon, Nov 16, 2009 at 11:25 AM, MS Al-Rawi <raw...@yahoo.com> wrote:

Hi ....... In a previous post, I noted that too many acts have NaN values when running GNB. I was reviewing the GNB code and I noticed that the joint probability is used without normalization. I quote the following comment from the file 'test_gnb.m'  % with uniform priors, the final log posterior is just the  % normalized log likelihood.  We would normally try to normalize  % within logs to avoid underflow, but that doesn't seem to be  % necessary here.  This is true when very few voxels are used, but when thousands  (and may be hundreds) of voxels  are used as input, it will lead to zeroing the posteriori probability. I think it is better (and might be necessary)  to do the normalization as follows (to be placed right after the above comment):  ZZ = abs(sum(log_likelihood,1)/nConds); % I found this to be very useful  for k=1:nTimepoints      log_likelihood(:,k) = log_likelihood(:,k)/ZZ(k);  end  I think, however, that the above is not 100% correct since the normalization is placed before exp(log_likelihood) which violates Bayes formula. Thus, to solve the problem we may just use the log_likelihood itself, after adding a simple shifting statement. I copy paste 'test_gnb' after slightly modifying the code to the one I propose. Regards Al-Rawi.. function [acts,scratch] = test_gnb(testpats,testtargs, scratch) % Use a Gaussian Naive Bayes classifier to predict regressors. % % [ACTS, SCRATCHPAD] = TEST_GNB(TESTPATS, TESTTARGS, SCRATCH) % % See TRAIN_GNB for more info. % License: %===================================================================== % % This is part of the Princeton MVPA toolbox, released under % the GPL. See http://www.csbmb.princeton.edu/mvpa for more % information. % % The Princeton MVPA toolbox is available free and % unsupported to those who might find it useful. We do not % take any responsibility whatsoever for any problems that % you have related to the use of the MVPA toolbox. % % ====================================================================== nConds = size(testtargs,1); [nVox nTimepoints] = size(testpats); % To make a prediction for a given test pattern, we compute the % posterior likelihood of observing the label for the given % pattern: % % Pr(Y = y | X = x) ~ Pr( X = x | Y = y) * Pr (Y = y) warning('off'); % compute the likelihood of the data under the MLE estimated % gaussian model for each category for k = 1:nConds % we calculate the proper mu and sigma for each voxel and each % timepoint: % scratch.mu is a nVox x 1 vector of voxel means for that condition mu = repmat(scratch.mu(:,k), 1, nTimepoints); % scratch.sigma is a nVox x 1 vector of voxel variances for a condition sigma = repmat(scratch.sigma(:,k), 1, nTimepoints); % compute the likelihood of the data under label k (this just the % value of the gaussian equation using the estimated means and variances) raw_likelihood = normpdf(testpats, mu, sigma); % GNB assumes all voxels are independent, so we multiply these probabilities % together to get the likelihood of each data point: we do this in % the log domain by summing to avoid underflow log_likelihood(k, :) = sum(log(raw_likelihood), 1); end % makes it more correct if we want to add a prior log_posterior = log_likelihood + repmat(log(scratch.prior), 1, nTimepoints); % The exp will result in underflow as in acts = exp(log_posterior), we might face cases like exp(-1000) which gives 0. However, % if P(w1) > P(w2), then log(P(w1)) > log(P(w2)) , so lets forget the exp and normalization and use the log itself posterior = log_posterior - repmat( min(log_posterior),nConds,1) ; % simple shifting to produce +ve values acts = posterior; acts = acts ./ repmat(sum(acts,1), nConds, 1);

[MS Al-Rawi]

Thank you Francisco for your response, there was a problem in the layout, thus, I had to post and remove few times. Anyway, as for what you mentioned

"but what you describe sounds like what people generally do: they work in log-space, without using the normalization constant."

Thats quite right, however, the underflow problem occurs when we try to use the exp() function to restore the original posteriori probability value which. I suggest to remove this function call, and I hope that this removal causes no problem. We may use the following form which performs simple shifting: 

posterior = log_posterior - repmat( min(log_posterior),nConds,1) ; % simple shifting to produce +ve values

Example: Using the subject provided by MVPA (Haxby 2000 et al), there were 520 underflow values leading to NaN (not a number) out of (8x121 =968) posteriori's (this subject has only ~577 masked voxels and the problem is worse if more voxels are used). One cannot notice this unless he debugs into the code line by line, and its not enough to judge by looking at the final classification result(s).

If any concern regarding the above shifting and exp removal let me know please so that we find another treatment.

All the best

Rawi

---

From: Francisco Pereira <francisc...@gmail.com>
To: mvpa-t...@googlegroups.com
Sent: Mon, November 16, 2009 4:32:14 PM
Subject: [mvpa-toolbox] Re: Revisiting the issue of Gaussian Naive Bayes classifier- new and Final one.... I hope!

--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups "Princeton MVPA Toolbox for Matlab" group.
To post to this group, send email to mvpa-t...@googlegroups.com
To unsubscribe from this group, send email to mvpa-toolbox...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/mvpa-toolbox?hl=en
-~----------~----~----~----~------~----~------~--~---