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. We have some sense that a small ${\chi }^2$ means a good fit and a big one means a bad fit, but how small and how big? How good and how bad? A ${\chi }^2$ statistic must be interpreted in the context of the relative sizes of $N$ and $M$. Obviously, if you allow $M=N$, then you have as many degrees of freedom in your model as there are observations to test its fit, and so it's trivial to get ${\chi }^2=0$ in this case.

Imagine that we knew the true, measurement-error-free value of every observation; let's call that ${\overline{y}}_i$. Imagine also that our reckoning of the measurement error distributions (that is, that each one is Gaussian with mean ${\overline{y}}_i$ and standard deviation ${\sigma }_i$) is exactly correct. Then each noisy observation $y_i$ is made from a Gaussian distribution, and we can normalize this observation to a z-score (standard normal variate) by $z_i=(y_i-{\overline{y}}_i)/{\sigma }_i$. Note, then, that the sum of the squares of all the $z_i$ is a ${\chi }^2$ statistic. The sum of the squares of $k$ independent standard normal random variables itself forms a probability distribution, and this distribution has a name: the ${\chi }^2$ distribution with $k$ degrees of freedom.

Now flip this concept around: rather than knowing true, error-free ${\overline{y}}_i$ values and talking about the probability distribution of the $y_i$ values, what we know is the $y_i$s and we want to evaluate the probability that our model has estimated the ${\overline{y}}_i$s correctly. The sum of the squares of the normalized differences between $y(x_i|a)$ and $y_i$ should follow the same ${\chi }^2$ distribution as above, if the model estimated ${\overline{y}}_i$ correctly. For every degree of freedom in our model ($M$), though, we have to eliminate one degree of freedom (that is, consideration of one observation) from the ${\chi }^2$ distribution, since theoretically we could've devoted that one degree of freedom to exactly reproducing one noisy observation that we made.

What we've arrived at is ${\chi }^2$ model selection. Given a value ${\chi }^2$ as computed above, the probability $Q$ of seeing an error as bad or worse than the observed ${\chi }^2$ under a different set of observations drawn from the same underlying distribution, given the degrees of freedom $\nu \equiv N-M$, is:

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

where $f(\psi |\nu )$ is the PDF of the ${\chi }^2$ distribution and $F(\psi |\nu )$ is of course its CDF. 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 it's almost impossible to get a better fit to the observations, even given the number of degrees of freedom. High values are actually sort of suspicious, and either indicate that the ${\sigma }_i$s were overestimated, $M$ was underestimated, or the data were actually fudged to get a better result.

This means that maximizing $Q$ over all values of $M$ (from very simple to very complex models) should find the "happy medium" model that maximizes explanatory power (low ${\chi }^2$) while minimizing model complexity (low $M$).

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

Bundle Adjustment

• $N$ is the number of voxels in all the DWIs that should be included in the model. A given voxel should be included if it is part of a liberal white-matter mask or has any part of the bundle model passing through it.
  • The point of this formalism is to avoid including background voxels unnecessarily; we could theoretically drown out the effects of anything interesting by creating an arbitrarily large FOV that included a huge volume of empty space. Limiting the definition of voxels that we want to include in the goodness-of-fit measure keeps us honest and sensitive.
• $M$ is the total number of degrees of freedom in the bundle model, which are difficult to count. They include at least: 3 DOF per triangulation vertex, 3 DOF per triangle (vertex indices), 1 DOF per bundle (the newline separating lists of triangles from each other), and some DOF for the internal "microstructure" model in each bundle.
• ${\sigma }_i$ is the background noise level. We use it (see the section on microstructure fitting) to convert Rician-distributed random variables to equivalent standard Gaussian RVs.

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 noise in the background of 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}}}_{\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{e}}(x|y)\mathrm{d}x=\underset{-\mathrm{\infty }}{\overset{{\tilde{y}}_i}{\int }}{\mathbb{\mathbb{P}}}_{\mathrm{g}\mathrm{a}\mathrm{u}\mathrm{s}\mathrm{s}}(x|y)\mathrm{d}x$

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.