pyhctsa.operations.spectral.sinusoid_fit

pyhctsa.operations.spectral.sinusoid_fit(y, model='sin1')

Fit sinusoids or a Fourier series to the time series.

Fits a sum of 1-3 sinusoids, or a Fourier series with 1-3 terms, to the time series as a function of its time index t = 0..N-1. The values are fitted in the order in which they occur, so the result depends on the temporal ordering of the data. (This is the time-series-model branch of the former DN_SimpleFit, split off in hctsa as SP_SinusoidFit because the distribution of values is unaffected by temporal ordering, whereas these fits are.)

The fitted models are:

  • 'sinK': a sum of K sinusoids, sum_i a_i*sin(2*pi*f_i*t + c_i), with free amplitudes, phases, and frequencies f_i in [1/N, 1/2 - 1/(2N)] cycles per sample. The amplitudes and phases are linear parameters, found by least squares for given frequencies, and the frequencies are searched deterministically (no random starts: a grid search refined by Gauss-Newton iterations): pyhctsa.robust.bf_fit_sinusoids(). The frequencies are bounded below by 1/N, so that each sinusoid completes at least one full cycle within the series: these models measure how well the series is described by oscillations, and a component with less than one cycle in the window is a trend (its period is not observed), which would let a sinusoid fit a trending series as if it were periodic (and below 1/(2N) its amplitude and phase are not determined: it cannot be told apart from a constant plus a linear trend).

  • 'fourierK': a K-term Fourier series, a0 + sum_i (a_i*cos(i*w*t) + b_i*sin(i*w*t)), with a single fitted fundamental frequency w, so that the terms are harmonically related.

Goodness of fit is summarized by the root mean square error (and R^2), and the residuals are characterized by their autocorrelation at lags 1 and 2 and by a runs test, which together reveal remaining temporal structure that the model has not captured (e.g., whether a periodic component has been fully explained).

Parameters:
y : array-like

The input time series.

model : {'sin1', 'sin2', 'sin3', 'fourier1', 'fourier2', 'fourier3'}, optional

The model to fit. Default is 'sin1'.

Returns:

  • r2: the R^2 of the fit, relative to the mean of the series. The 'sinK' models have no intercept (constant) term, so r2 can be negative for them.

  • adjr2: the degrees-of-freedom-adjusted R^2.

  • rmse: the root mean square error of the fit (the residual sum of squares divided by the degrees of freedom of the error).

  • resAC1, resAC2: the autocorrelation of the residuals at lags 1 and 2 ('Fourier' method).

  • resrunsz: the signed z-statistic of a runs test on the residuals (pyhctsa.robust.bf_runs_z()): negative when the residuals have fewer runs about their median than expected for a random order (slowly varying residuals).

The three residual outputs are NaN if the residuals have no variation beyond floating-point round-off, as for an exact fit (pyhctsa.robust.bf_residual_stats()). NaN is returned instead of a dict if the model cannot be fitted (or, for K sinusoids, if there are fewer than 3K + 1 samples).

Return type:

dict or float

Notes

The Fourier series are fitted by nonlinear least squares (the linear coefficients and the fundamental frequency) with SciPy’s trust-region-reflective least_squares, started from the periodogram peak of the data (see _fourier_start_point()). The fit is multimodal in the frequency, so the result depends on the start point.