User:Jadrian Miles/Thesis manifesto: probabilistic worldview

My dissertation involves solving three problems:

1. Automatically clustering tractography curves together so that the resulting clusters are neither too small (low reconstruction error, high model complexity) nor too big (high reconstruction error, low model complexity)
2. Automatically adjusting macrostructure elements to match input DWIs so that the elements' surfaces are neither too bumpy (low reconstruction error, high model complexity) nor too smooth (high reconstruction error, low model complexity)
3. Automatically adjusting microstructure properties within a given region in space to match input DWIs so that the spatial frequency of the microstructure parameters is neither too high (low reconstruction error, high model complexity) nor too low (high reconstruction error, low model complexity)

In each case, our goal is to balance the tradeoff between reconstruction error and model complexity. In order to do so in a principled fashion, we must define, for each case, what is meant by reconstruction error and model complexity. What probabilistic tools are available for us to go about this?

${\chi }^2$ Fitting

One apparent option is ${\chi }^2$ fitting. Given:

• $N$ observations $y_i$ at independent variables $x_i$
• The knowledge that the measurement error on each observation is Gaussian with standard deviation ${\sigma }_i$
• A model instance with $M$ parameters $a_j$, which gives a reconstructed "observation" for a given independent variable as $y(x_i|a_0\dots a_{M-1})$

we define the ${\chi }^2$ statistic as follows:

${\chi }^2=\sum _{i=0}^{N-1}{\left(\frac{y_i-y(x_i|a_0\dots a_{M-1})}{{\sigma }_i}\right)}^2$

This is essentially a normalized error score, and in itself tells us nothing about goodness of fit. The interesting thing that we can do here, though, is compute the probability of seeing an error as bad or worse than the observed error score under a different set of observations drawn from the same underlying distribution, given the degrees of freedom $\nu \equiv N-M$. This probability is called $Q$, and is defined as:

$Q=1-\underset{{\chi }^2}{\overset{\mathrm{\infty }}{\int }}F(\psi ,\nu )\mathrm{d}\psi$

where $F(\psi ,\nu )$ is the CDF of the ${\chi }^2$ distribution. Here's how to interpret $Q$:

• Really low ($Q<0.001$) means that you can't do a much worse job describing the data than you did with your model instance. This means either that your fit is bad (${\chi }^2$ is really high) or that your model has so many degrees of freedom (and thus $\nu$ is so small) that even a low ${\chi }^2$ isn't convincing.
• Really high ($Q\approx 1$) means that you can't do much better, even given the number of degrees of freedom. This is a suspicious result and may indicate data-fudging or just overestimation of the ${\sigma }_i$s, so it shouldn't come up in an automated system.

This means that maximizing $Q$ over all values of $M$ (from very simple to very complex models) should find the "happy medium".

Let's consider the meaning of all these variables for each of the problems defined above:

Curve Clustering

• $N$ is the number of vertices in all the curves in the tractography set. Each $x_i$ therefore represents the selection of a single vertex (by two parameters: the index of the curve in the curve set, and the index or arc-length distance of the desired vertex along this curve). $y_i$ isn't particularly well-defined in this case (see below), but amounts to the description of the curve near $x_i$: the vertex's position and the angle between the consecutive segments, maybe.
• $M$ is the total number of curves that form the skeletons of the shrink-wrap polyhedra for the model, plus the total number of bundles/polyhedra. Note that even in the most "complicated" model, $M$ is smaller than $N$ by a factor of half the average number of $x_i$s per curve. $y(x_i|a_0\dots a_{M-1})$ is the description of the reconstructed curve specified by $x_i$ at the specified position along it, given the association of curves with bundles specified by the model parameter values $a_0\dots a_{M-1}$.
• ${\sigma }_i$ is defined according to the definition of the difference $(y_i-y(x_i|a_0\dots a_{M-1}))$ that we choose (see below).

Let's assume that for comparing the reconstruction ($y(x_i|a_0\dots a_{M-1})$, which we will just call $y$) with the observation ($y_i$), we want to incorporate both a position and angle reconstruction error. Then we must define some sort of scalar pseudo-difference function

$f(y,y_i)=\alpha \times d_p(y,y_i)+(1-\alpha )\times d_{\theta }(y,y_i)$

where $d_p$ is the difference in positions (Euclidean distance) between the vertices, and $d_{\theta }$ is some difference between angles, perhaps one minus the dot product. $\alpha$ may either be hand-tuned or solved for in the course of the optimization; in this latter case, $M$ must be increased by one.

The ${\sigma }_i$s must be determined experimentally (with synthetic phantoms, for example). It remains for these experiments to demonstrate that the errors are normally distributed. If not, some wacky transform might be necessary to shoehorn it into the ${\chi }^2$ regime, or Monte-Carlo simulation of the error distribution might allow for a different but equivalent approach (see §15.6 of NR).

Microstructure Fitting

• $N$ is the number of voxels in the relevant region multiplied by the number of observations (diffusion weightings) in each. Each individual configuration (voxel and diffusion weighting) is an $x_i$. The observed intensity for each $x_i$ is $y_i$.
• $M$ is the number of control parameters in the space-filling microstructure model. Assuming that the microstructure instance for each voxel is determined by some sort of spline with sparse control points, $M$ is 3 parameters for position plus the number of parameters actually controlling the microstructure model, multiplied by the number of control points. $y(x_i|a_0\dots a_{M-1})$ is the reconstructed signal intensity in the voxel and diffusion weighting specified by $x_i$, given the model parameter values $a_0\dots a_{M-1}$.
• ${\sigma }_i$ is known for all $i$ by estimation of the Rician noise parameters from the DWIs. But...

Unfortunately, we know that the noise on the $y_i$s is not Gaussian but rather Rician. Therefore $(y_i-y(x_i|a))$ is not normally-distributed, and so the ${\chi }^2$ statistic is invalid in the form above. We can, however, create a "fake" observation ${\tilde{y}}_i$ from $y_i$ for which $({\tilde{y}}_i-y(x_i|a))$ is normally distributed. We do this by solving for the ${\tilde{y}}_i$ that satisfies the following equality:

$\underset{-\mathrm{\infty }}{\overset{y_i}{\int }}\mathbb{\mathbb{P}}(y_i|y)\mathrm{d}{y}_i=\underset{-\mathrm{\infty }}{\overset{{\tilde{y}}_i}{\int }}\mathbb{\mathbb{P}}({\tilde{y}}_i|y)\mathrm{d}{\tilde{y}}_i$

The two sides of the equality are values of the cumulative distribution functions of, respectively, a Rician random variable and a gaussian random variable:

$1-Q_1(\frac{y}{{\sigma }_i},\frac{y_i}{{\sigma }_i})=\frac{1}{2}(1+\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{{\tilde{y}}_i-y}{\sqrt{2{\sigma }_i^2}}\right))$

where $Q_1$ is the first-order Marcum Q-function, and $\mathrm{e}\mathrm{r}\mathrm{f}$ is the Gauss error function. Solving, we find:

${\tilde{y}}_i=y+\sqrt{2{\sigma }_i^2}{\mathrm{e}\mathrm{r}\mathrm{f}}^{-1}(1-2Q_1(\frac{y}{{\sigma }_i},\frac{y_i}{{\sigma }_i}))$

This function re-maps a Rician random variable $y_i$ into a Gaussian random variable ${\tilde{y}}_i$ with the same ${\sigma }_i$ and central value $y$, and so we may use this ${\tilde{y}}_i$ in place of $y_i$ in the summation for the ${\chi }^2$ statistic.

The inner loop of our algorithm is given a model with $M$ degrees of freedom, and fits the $a_j$s to maximize $Q$ (or equivalently, minimize ${\chi }^2$) for that model (see §15.5 of NR). The outer loop is a search over all models (some of which may have the same $M$) to find a global maximum in $Q$. As mentioned above, we should never get an unrealistically high $Q$ value, as that's really only achievable by either fudging the data or overestimating the ${\sigma }_i$s.