Implementation

Both forms of priors lead to posteriors that depend on $\beta$ and possibly $\lambda$ where these cannot be easily integrated over analytically. Therefore the alternatives are: (1) to approximate the integration numerically; (2) to simplify the models/assumptions (e.g. flat priors on $\alpha$); or (3) to set $\beta$ and $\lambda$ to be known constants (pragmatically they can be measured from the data $Y$).

The most expensive computation is that of the residuals, as this requires the accumulation of intensities over many voxels and the appropriate updating of summary statistics to do the effective planar fit. Once these statistics have been generated, the remaining matrix computations are relatively fast and so it is feasible to integrate over $\beta$ numerically, as the residual term is easily and cheaply recalculated.