User:Jadrian Miles/Thesis manifesto: probabilistic worldview

boopboopboopMy dissertation involves solving three problems: boopboopboop# 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) boopboopboop# 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) boopboopboop# 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) boopboopboop boopboopboopIn 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? boopboopboop boopboopboop==${\chi }^2$ Fitting== boopboopboop boopboopboopOne apparent option is ${\chi }^2$ fitting. Given: boopboopboop* $N$ observations $y_i$ at independent variables $x_i$ boopboopboop* The knowledge that the measurement error on each observation is Gaussian with standard deviation ${\sigma }_i$ boopboopboop* 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})$ boopboopboopwe define the ${\chi }^2$ statistic as follows: boopboopboop boopboopboop${\chi }^2=\sum _{i=0}^{N-1}{\left(\frac{y_i-y(x_i|a_0\dots a_{M-1})}{{\sigma }_i}\right)}^2$ boopboopboop boopboopboopThis 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. boopboopboop boopboopboopImagine 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. boopboopboop boopboopboopNow 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. boopboopboop boopboopboopWhat 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: boopboopboop boopboopboop$Q({\chi }^2|\nu )=\underset{{\chi }^2}{\overset{\mathrm{\infty }}{\int }}f(\psi ,\nu )\mathrm{d}\psi =1-F({\chi }^2|\nu )$ boopboopboop boopboopboopwhere $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$: boopboopboop* 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. boopboopboop* 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. boopboopboop boopboopboopThis 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$). boopboopboop boopboopboopLet's consider the meaning of all these variables for each of the problems defined above: boopboopboop boopboopboop===Curve Clustering=== boopboopboop* $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. boopboopboop* $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}$. boopboopboop* ${\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). boopboopboop boopboopboopLet'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 boopboopboop boopboopboop$f(y,y_i)=\alpha \times d_p(y,y_i)+(1-\alpha )\times d_{\theta }(y,y_i)$ boopboopboop boopboopboopwhere $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. boopboopboop boopboopboopThe ${\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). boopboopboop boopboopboop===Bundle Adjustment=== boopboopboop boopboopboop boopboopboop===Microstructure Fitting=== boopboopboop* $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$. boopboopboop* $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}$. boopboopboop* ${\sigma }_i$ is known for all $i$ by estimation of the Rician noise parameters from the DWIs. But... boopboopboop boopboopboopUnfortunately, 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: boopboopboop boopboopboop$\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$ boopboopboop boopboopboopThe two sides of the equality are values of the cumulative distribution functions of, respectively, a Rician random variable and a Gaussian random variable: boopboopboop boopboopboop$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))$ boopboopboop boopboopboopwhere $Q_1$ is the first-order Marcum Q-function, and $\mathrm{e}\mathrm{r}\mathrm{f}$ is the Gauss error function. Solving, we find: boopboopboop boopboopboop${\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}))$ boopboopboop boopboopboopThis 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. boopboopboop boopboopboopThe 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. boopboopboop boopboopboop= Notes =

boopboopboop