ThomasYeoLab / CBIG

MIT License
570 stars 380 forks source link

Potential algorithm implementation error on group parcellation in MS-HBM (a.k.a. Yeo 2011 parcellation method) #56

Closed shuqike closed 1 month ago

shuqike commented 1 month ago

Expected situation

According to the paper "Clustering on the Unit Hypersphere using von Mises-Fisher Distributions" by Arindam Banerjee et al., the normalizing constant of a von Mises-Fisher distribution is given by

$$ cd(\kappa)=\frac{\kappa^{d/2-1}}{(2\pi)^{d/2}I{d/2-1}(\kappa)} $$

where $d$ is the dimension of a feature vector.

Actual situation

However, in the code CBIG_VonmisesSeriesClustering_fix_bessel_randnum_bsxfun.m line 129: tic, clustered = direcClus_fix_bessel_bsxfun(series, num_clusters, size(series, 2) - 1, num_tries, lambda, 0, 0, 1e-4, 1, max_iter, 1); you passed size(series, 2) - 1 into the direcClus_fix_bessel_bsxfun function. Is this equivalent to $c_{d-1}(\kappa)$?

In this case, should we use size(series, 2) without minus 1 instead?

Thanks!

rubykong commented 1 month ago

@shuqike you are right that actually no need to -1. This is a very old bug thomas made long time ago and we didn't want to fix. Because this only affects the approximation. I am referring to AppendixB in Lashkari et al., NeuroImage 2010 (Discovering structure in the space of fMRI selectivity profiles). We are using the upper bound as the final approximation. However, using D-1 rather than D still holds the approximation (Equation B.7). It's just that it's no longer the upper bound. We don't want to fix it because it will affect the replication of some old scripts. And since it doesn't really affect the performance, we decide to keep it to D-1.

shuqike commented 1 month ago

I see. Thanks for quick response!