Skip to content

Getting started

This page fits a damped oscillator from noisy measurements: building the data, writing the equation, reading the result, and checking how sure the estimates are.

1. The data

pysnape takes measurements on a grid as a Data object: the fields you measured, the coordinates of the grid, and any known inputs such as an applied force. Here we make a damped oscillation, x = exp(-0.1 t) cos(2 t), and add 5 % noise:

import numpy as np
import pysnape as ps

t = np.linspace(0.0, 10.0, 301)
clean = ps.Data(fields={"x": np.exp(-0.1 * t) * np.cos(2.0 * t)}, coords={"t": t})
data = clean.add_noise(percent=5, seed=0)

With your own measurements, put your arrays in fields and your time (or space) samples in coords. A field's array has one axis per coordinate, in the order the coordinates are given. Data covers files, missing points, measured accelerations and the published datasets.

2. The equation

Write it as on paper. Derivatives are suffixes (x_t, x_tt, u_xx), and every name that is not a field, a coordinate, a known input or a known constant is an unknown coefficient, which needs a starting value:

equation = "x_tt + a*x_t + b*x = 0"
start = {"a": 1.0, "b": 1.0}

The starting values need not be close: by default pysnape also tries a grid of starting points scaled to the data. Writing equations lists everything an equation may contain.

3. Fit

result = ps.fit(equation, data, unknowns=start)
print(result.summary())

The summary lists the equation, the basis it chose, each coefficient's start and estimate, and three numbers that say how the fit went:

  • data misfit: how far the fitted signal is from the measurements, relative to their size. For noisy data it is about the noise level (here about 0.05); for clean data it should be small.
  • equation residual: how far the fitted signal is from satisfying the equation.
  • equation weight: how strongly the equation was enforced in the answer chosen (How the solver works).

The truth here is a = 0.2, b = 4.01:

print(result.theta)

4. Use the result

x_fit = result.field()  # the fitted signal on the data grid
v_fit = result.field(deriv="t")  # its time derivative, no differencing of noisy data
r = result.residual()  # the equation residual on the grid

result.plot() draws the measurements against the fitted signal (it needs matplotlib: pip install "pysnape[plot]").

5. How sure is it?

For measured data, refit with resampled residuals (a residual bootstrap):

boot = ps.bootstrap(result, n=20, jobs=2)
print(boot.summary())  # mean, coefficient of variation and 95 % interval per coefficient

For a synthetic study, where you have the clean signal, refit fresh noise draws instead:

truth_fit = ps.fit(equation, clean, unknowns=start, multistart=False)
ens = ps.ensemble(truth_fit, n=20, noise=5, jobs=2)
print(ens.summary(truth={"a": 0.2, "b": 4.01}))

Uncertainty explains the difference.

When it goes wrong

  • A warning that the basis is probably too coarse. The fitted signal misses the data by more than 25 %. If your noise is not that large, the spline cannot follow the solution under the equation: give knots= more knots in the coordinate where the solution changes fastest, or a list of candidates, knots=[100, 200, 400], and pysnape keeps the best.
  • A slow fit. Pass verbose=True to see each stage as it finishes, and see the speed notes in How the solver works.
  • An error about the equation. The message names the term and the reason; the equation rules list what is allowed.