9. Cramér-Rao Bounds

DRSuite provides a tool to evaluate estimation-theoretic Cramér-Rao bounds for this type of spectrum estimation problem. This provides a quantitative measure that can be used for evaluating the quality of estimated spectra, and can also be used to prospectively assess the estimation-theoretic characteristics of different experimental acquisition paradigms. This can be very useful for designing and optimizing experiments, as described in the following references (which should be cited, along with DRSuite, when using this tool):

    1. Kim, J. L. Wisnowski, C. T. Nguyen, J. P. Haldar. Multidimensional Correlation Spectroscopic Imaging of Exponential Decays: From Theoretical Principles to In Vivo Human Applications. NMR in Biomedicine 33:e4244, 2020. <https://doi.org/10.1002/nbm.4244>

    1. Kim, J. L. Wisnowski, J. P. Haldar. Improved Efficiency for Microstructure Imaging using High-Dimensional MR Correlation Spectroscopic Imaging. Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, 2017, pp. 1264-1268. <https://doi.org/10.1109/ACSSC.2017.8335555>

    1. Kim, J. P. Haldar. Faster Diffusion-Relaxation Correlation Spectroscopic Imaging (DR-CSI) using Optimized Experiment Design. International Society for Magnetic Resonance in Medicine 25th Annual Meeting, Honolulu, 2017, p. 176. <https://archive.ismrm.org/2017/0176.html>

The software currently calculates Cramér–Rao bounds for a single-voxel multi-compartment signal model. Given:

  • a parametric multicomponent signal expression (e.g., \(d(b,TE) = \sum_{s=1}^S f_s \exp(-b D_s) \exp(-TE/T2_s)\));

  • a list of spectral parameters (e.g., \(D\), \(T2\));

  • a list of acquisition parameters that influence image contrast (e.g., \(b\), \(TE\));

  • the acquisition parameters sampled by a given experiment (e.g., \(\{(b_i,TE_i)\}_{i=1}^{N_a}\));

  • the parameters corresponding to the desired “true” model to evaluate

The tool builds the Fisher information matrix for this model and computes the Cramér-Rao bound by taking the diagonal of the inverse of the Fisher information matrix. The computation can be performed with the following command:

estimate_crlb.sh -funcfile func_expression.txt -outprefix Result/CRLB_values

estimate_crlb('funcfile','func_expression.txt', ...
          'outprefix','Result/CRLB_values')

9.1. Required Inputs

  • funcfile

    Filename of a .txt file storing the signal model and its associated parameters required for Cramér-Rao bbound calculation. The file must define the following variables (with these exact names):

    func        = '';
spect_params = '';
acq_params   = '';
exp_vals     = [];
components   = [];

    The meaning of each field is:

    • func A MATLAB-style function expression for the noise-free signal for a single component, written in terms of the spectral parameters and acquisition parameters.

      Example (diffusion–T2 model):

      func = '(1-exp(-b*D))*(exp(-TE/T2))';

    • spect_params Comma-separated list of spectral parameter names that appear in func. These names must match the symbols used in func. Example:

      spect_params = 'D,T2';

    • acq_params Comma-separated list of acquisition parameter names that appear in func. The order of this list must match the column order of exp_vals.

      Example:

      acq_params = 'b,TE';

    • exp_vals A matrix, where each row stores one combination of acquisition parameters to be used when building the Fisher information matrix. Columns must follow the order specified in acq_params.

      Example:

      exp_vals = [ ...
    0,   99;
    0,  120;
    0,  160;
    0,  200;
    0,  250;
    0,  300;
    0,  400;
    0.2, 99;
    ...
    5,  400 ];

      In this example, the first column corresponds to b and the second column to TE.

    • components A matrix of size \(S \times (M+1)\), where each row specifies the prior information of one spectral component (total \(S\) components). The columns must be \([f_s, \text{spect_param1}_s, \text{spect_param2}_s, \cdots, \text{spect_paramM}_s]\). The spectrum parameters in each row, except \(f_s\), must follow the same order specified in spect_params. Example for two components:

      components = [1, 0.5, 70;    % f1, D1, T2_1
              1, 1.3, 100];  % f2, D2, T2_2

  • outprefix

    String specifying the filename prefix where the CRLB results will be stored. The tool appends .txt automatically. For example, outprefix = 'Result/CRLB_values' will create a file named Result/CRLB_values.txt.

9.2. Optional Field Inside funcfile

  • noise_std (scalar, optional)

    Standard deviation of the additive white Gaussian noise assumed in the model. If not provided, the tool uses

    noise_std = 1;

    The Fisher information scales as \(1/\text{noise_std}^2\) and therefore the CRLB values scale as \(\text{noise_std}^2\). Use a realistic noise_std consistent with the scale of your signal.

9.3. Output

The tool writes a single tab-separated text file:

  • File: [outprefix].txt

  • Size: \(S \times (M+1)\).

Each row corresponds to one spectral component, and each column corresponds to the variance lower bound (CRLB) of one parameter:

\[\text{Row } s = \big[ \mathrm{CRLB}(w_s),\ \mathrm{CRLB}(\text{param1}_s),\ \mathrm{CRLB}(\text{param2}_s) \big].\]

Notes:

  • The CRLB values have the units of variances:

    • \(\mathrm{CRLB}(w_s)\) → (units of weight)²

    • \(\mathrm{CRLB}(D_s)\) → (units of \(D\))²

    • \(\mathrm{CRLB}(T_{2,s})\) → (units of \(T_2\))²

  • To obtain the minimum achievable standard deviation for each parameter, take the square root of the corresponding CRLB entry.

  • For physically meaningful results, ensure that the units used in func, exp_vals and components are consistent (e.g., TE and T2 both in ms, b and D chosen so that b*D is dimensionless).