Struct varpro::statistics::FitStatistics

pub struct FitStatistics<Model>
where
    Model: SeparableNonlinearModel,
    DefaultAllocator: Allocator<Model::ScalarType, Dyn, Dyn> + Allocator<Model::ScalarType, Dyn>,
{ /* private fields */ }

This structure contains some additional statistical information about the fit, such as errors on the parameters and other useful information to assess the quality of the fit.

Where is $R^2$?

We don’t calculate $R^2$ because “it is an inadequate measure for the goodness of the fit in nonlinear models” (Spiess and Neumeyer 2010). See also here, here, and here, where the recommendation is to use the standard error of the regression instead. See also the next section.

Model Selection

If you want to want to use goodness of fit metrics to decide which of two models to use, please look into the Akaike information criterion, or the Bayesian information criterion. Both of these can be calculated from standard error of the regression. Note that this fit uses least squares minimization so it assumes the errors are zero mean and normally distributed.

Implementations

impl<Model> FitStatistics<Model>

where Model: SeparableNonlinearModel, DefaultAllocator: Allocator<Model::ScalarType, Dyn, Dyn> + Allocator<Model::ScalarType, Dyn>,

pub fn covariance_matrix(&self) -> &OMatrix<Model::ScalarType, Dyn, Dyn>

The covariance matrix of the parameter estimates. Here the parameters are both the linear as well as the nonlinear parameters of the model. The linear parameters are ordered first, followed by the non-linear parameters as if we had one combined parameter vector $(\vec{c}, \vec{\alpha})$. See O’Leary and Rust 2012 for reference.

Example

Say our model has two linear coefficients $\vec{c}=(c_1,c_2)^T$ and three nonlinear parameters $\vec{\alpha}=(\alpha_1,\alpha_2,\alpha_3)^T$. Then the covariance matrix is odered for the parameter vector $(c_1,c_2,\alpha_1,\alpha_2,\alpha_3)^T$. The covariance matrix $C$ (upper case C) is a square matrix of size $5 \times 5$. Matrix element $C_{ij}$ is the covariance between the parameters at indices $i$ and $j$, so in this example:

• $C_{11}$ is the variance of $c_1$,
• $C_{12}$ is the covariance between $c_1$and $c_2$,
• $C_{13}$ is the covariance between $c_1$ and $\alpha_1$,
• and so on.

pub fn correlation_matrix(&self) -> OMatrix<Model::ScalarType, Dyn, Dyn>

where Model::ScalarType: Float,

👎Deprecated since 0.9.0: Use the method calc_correlation_matrix.

calculate the correlation matrix. Deprecated, use the `calc_correlation_matrix`` function instead.

pub fn calculate_correlation_matrix( &self ) -> OMatrix<Model::ScalarType, Dyn, Dyn>

where Model::ScalarType: Float,

The correlation matrix, ordered the same way as the covariance matrix.

Note The correlation matrix is calculated on the fly from the covariance matrix when this function is called. It is suggested to store this matrix somewhere to avoid having to recalculate it.

pub fn weighted_residuals(&self) -> OVector<Model::ScalarType, Dyn>

the weighted residuals

\vec{r_w} = W * (\vec{y} - \vec{f}(\vec{\alpha},\vec{c}))

at the best fit parameters $\vec{alpha}$ and $\vec{c}$.

In case of a dataset with multiple members, the residuals are written one after the other into the vector

pub fn regression_standard_error(&self) -> Model::ScalarType

where Model::ScalarType: Float,

the regression standard error (also called weighted residual mean square, or sigma). Calculated as

\sigma = \frac{||\vec{r_w}||}{\sqrt{N_{data}-N_{params}-N_{basis}}},

where $N_{data}$ is the number of data points (observations), $N_{params}$ is the number of nonlinear parameters, and $N_{basis}$ is the number of linear parameters (i.e. the number of basis functions).

pub fn reduced_chi2(&self) -> Model::ScalarType

the reduced chi-squared after the fit, i.e. $\chi^2/\nu$, where $\nu$ is the number of degrees of freedom.

pub fn nonlinear_parameters_variance(&self) -> OVector<Model::ScalarType, Dyn>

where Model::ScalarType: Scalar + Zero, DefaultAllocator: Allocator<Model::ScalarType, Dyn>,

helper function to extract the estimated variance of the nonlinear model parameters. Those could also be manually extracted from the diagonal of the covariance matrix.

pub fn linear_coefficients_variance(&self) -> OVector<Model::ScalarType, Dyn>

where Model::ScalarType: Scalar + Zero, DefaultAllocator: Allocator<Model::ScalarType, Dyn>,

helper function to extract the estimated variance of the linear model parameters.

pub fn confidence_band_radius( &self, probability: Model::ScalarType ) -> OVector<Model::ScalarType, Dyn>

where Model::ScalarType: Float + One + Zero + CastF64,

Calculate the radius (i.e. half of the width) of the confidence band of the fitted model function at the best fit parameters evaluated at the support points.

This function was inspired by the function eval_uncertainty of the python package lmfit. It will produce an identical result (within numerical accuracy) for the same fit, given the same probability level.

Arguments

• probability the confidence level of the confidence interval, expressed as a probability value between 0 and 1. Note that it can’t be 0 or 1, but anything in between. Since least squares fitting assumes Gaussian error distributions, we can use the quantiles of the normal distribution to relate this to often used “number of sigmas”. For example the probability $68.3 \% = 0.683$ corresponds to (approximately) $1\sigma$.

Returns

The half-width of the confidence band with the given probability. To calculate the confidence band, we first need to calculate the fit result at the best fit parameters. Then the upper bound of the confidence interval is given by adding this radius and the lower bound is given by subtracting this radius.

Example

Fit model model to observations y and calculate the best fit and the 1-sigma confidence band upper and lower bounds.

use varpro::solvers::levmar::LevMarProblemBuilder;
use varpro::solvers::levmar::LevMarSolver;
let problem = LevMarProblemBuilder::new(model)
              .observations(y)
              .build()
              .unwrap();

let (fit_result,fit_statistics) = LevMarSolver::default()
              .fit_with_statistics(problem)
              .unwrap();

let best_fit = fit_result.best_fit().unwrap();
let cb_radius = fit_statistics.confidence_band_radius(0.683);
// upper and lower bounds of the confidence band
let cb_upper = best_fit + cb_radius;
let cb_lower = best_fit - cb_radius;

An Important Note on the Confidence Band Values

This library chooses to implement the confidence interval such that it gives the same results as the python library lmfit. That means that the confidence interval is given as if during the fitting process the weights had been uniformly scaled such that the reduced $\chi^2$ after fitting is equal to unity: $\chi_/nu^2 = \chi^2/\nu = 1$, where $\nu$ is the number of degrees of freedom (i.e. the number of data points minus the number of total fit parameters).

Scaling the weights does not influence the best fit parameters themselves, but it does influence the resulting uncertainties. Read here for an in-depth explanation. Briefly, the expected value for $\chi^2_\nu$ for a fit with $\nu$ degrees of freedom is one. Therefore it can be reasonable to apply the scaling such that we force $chi^2_\nu$ to take its expected value. This will correct an overestimation or underestimation of the errors and is often a reasonable choice to make.

However, this will make this confidence band inconsistent with the other information from the fit, such as the standard errors from the covariance matrix or the reduced $\chi^2$. Luckily, it’s easy to obtain the confidence bands with the errors exactly as given. Just multiply the result of this function with the reduced $\chi^2$ of this fit.

let (fit_result,fit_statistics) = LevMarSolver::default()
              .fit_with_statistics(problem)
              .unwrap();

let best_fit = fit_result.best_fit().unwrap();
// get the unscaled confidence band by multiplying with the
// reduced chi2
let cb_radius_unscaled = fit_statistics.reduced_chi2() * fit_statistics.confidence_band_radius(0.683);
// upper and lower bounds of the confidence band
let cb_upper = best_fit + cb_radius_unscaled;
let cb_lower = best_fit - cb_radius_unscaled;

Panics

Panics if probability is not within the half-open interval $(0,1)$.

Trait Implementations

impl<Model> Clone for FitStatistics<Model>

where Model: SeparableNonlinearModel + Clone, DefaultAllocator: Allocator<Model::ScalarType, Dyn, Dyn> + Allocator<Model::ScalarType, Dyn>, Model::ScalarType: Clone,

fn clone(&self) -> FitStatistics<Model>

Returns a copy of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<Model> Debug for FitStatistics<Model>

where Model: SeparableNonlinearModel + Debug, DefaultAllocator: Allocator<Model::ScalarType, Dyn, Dyn> + Allocator<Model::ScalarType, Dyn>, Model::ScalarType: Debug,

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.