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.

${\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. 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 you have so many degrees of freedom 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, 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:

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.

Let us consider a single observation in isolation. The observed signal magnitude is $y_i$ and the reconstructed value from the model is $y$. Assuming $y$ is a correct reconstruction, the marginal probability of observing $y_i$ is:

$\mathbb{\mathbb{P}}(y_i|y)\mathrm{d}{y}_i=\frac{1}{2\pi {\sigma }_i^2}\mathrm{e}\mathrm{x}\mathrm{p}(-(y^2+y_i^2)/2{\sigma }_i^2)\underset{-\pi }{\overset{\pi }{\int }}\mathrm{e}\mathrm{x}\mathrm{p}((y{y}_i/{\sigma }_i^2)\mathrm{c}\mathrm{o}\mathrm{s}\theta )\mathrm{d}\theta \mathrm{d}{y}_i$

(Clearly the same probability results if we instead reverse the roles of $y$ and $y_i$.)

We wish to define a normally-distributed "fake observation" ${\tilde{y}}_i$ with the same probability:

$\mathbb{\mathbb{P}}({\tilde{y}}_i|y)\mathrm{d}{\tilde{y}}_i=\frac{1}{\sqrt{2\pi {\sigma }_i^2}}\mathrm{e}\mathrm{x}\mathrm{p}(-({\tilde{y}}_i-y{)}^2/2{\sigma }_i^2)\mathrm{d}{\tilde{y}}_i=\mathbb{\mathbb{P}}(y_i|y)\mathrm{d}{y}_i$

Allowing $\mathrm{d}{\tilde{y}}_i=\mathrm{d}{y}_i$, we arrive after some algebraic manipulation at:

$({\tilde{y}}_i-y{)}^2=y^2+y_i^2-2{\sigma }_i^2(\underset{-\pi }{\overset{\pi }{\int }}\mathrm{e}\mathrm{x}\mathrm{p}((y{y}_i/{\sigma }_i^2)\mathrm{c}\mathrm{o}\mathrm{s}\theta )\mathrm{d}\theta -\frac{1}{2}\mathrm{l}\mathrm{o}\mathrm{g}(2\pi {\sigma }_i^2))$

This is just the numerator of the $i$-th term in the summation for the ${\chi }^2$ statistic comparing the model reconstruction $y$ to the (fake) observation $\tilde{y}$. The value of the integral depends on only one variable, $m\equiv y{y}_i/{\sigma }_i^2$, so it may be tabulated in a lookup table for fast computation.

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. 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.