02: Fitting Power Spectrum Models

Introduction to the module, beginning with the model object.

# Import the model object
from specparam import SpectralModel

# Import a utility to download and load example data
from specparam.utils.download import load_example_data

# Download example data files needed for this example
freqs = load_example_data('freqs.npy', folder='data')
spectrum = load_example_data('spectrum.npy', folder='data')

Model Object

At the core of the module is the SpectralModel object, which holds relevant data and settings as attributes, and contains methods to run the algorithm to parameterize neural power spectra.

The organization and use of the model object is similar to scikit-learn:

• A model object is initialized, with relevant settings

• The model is used to fit the data

• Results can be extracted from the object

Calculating Power Spectra

The SpectralModel object fits models to power spectra. The module itself does not compute power spectra. Computing power spectra needs to be done prior to using the specparam module.

The model is broadly agnostic to exactly how power spectra are computed. Common methods, such as Welch’s method, can be used to compute the spectrum.

If you need a module in Python that has functionality for computing power spectra, try NeuroDSP.

Note that model objects require frequency and power values passed in as inputs to be in linear spacing. Passing in non-linear spaced data (such logged values) may produce erroneous results.

Fitting an Example Power Spectrum

The following example demonstrates fitting a power spectrum model to a single power spectrum.

# Initialize a model object
fm = SpectralModel()

# Set the frequency range to fit the model
freq_range = [2, 40]

# Report: fit the model, print the resulting parameters, and plot the reconstruction
fm.report(freqs, spectrum, freq_range)

WARNING: Lower-bound peak width limit is < or ~= the frequency resolution: 0.50 <= 0.49
        Lower bounds below frequency-resolution have no effect (effective lower bound is the frequency resolution).
        Too low a limit may lead to overfitting noise as small bandwidth peaks.
        We recommend a lower bound of approximately 2x the frequency resolution.

==================================================================================================

                                       POWER SPECTRUM MODEL

                       The model was fit with the 'spectral fit' algorithm
               Model was fit to the 2-40 Hz frequency range with 0.49 Hz resolution

                               Aperiodic Parameters ('fixed' mode)
                                        (offset, exponent)
                                         -21.6185, 0.7160

                         Peak Parameters ('gaussian' mode) 3 peaks found
                                CF:   9.36, PW:   1.04, BW:   1.59
                                CF:  11.17, PW:   0.23, BW:   2.88
                                CF:  18.25, PW:   0.33, BW:   2.85

                                    Model fit quality metrics:
                                     error (mae) is 0.0356
                                  gof (rsquared) is 0.9829

==================================================================================================

Fitting Models with ‘Report’

The above method ‘report’, is a convenience method that calls a series of methods:

• fit(): fits the power spectrum model

• print_results(): prints out the results

• plot(): plots the data and model fit

Each of these methods can also be called individually.

# Alternatively, just fit the model with SpectralModel.fit() (without printing anything)
fm.fit(freqs, spectrum, freq_range)

# After fitting, plotting and parameter fitting can be called independently:
# fm.print_results()
# fm.plot()

WARNING: Lower-bound peak width limit is < or ~= the frequency resolution: 0.50 <= 0.49
        Lower bounds below frequency-resolution have no effect (effective lower bound is the frequency resolution).
        Too low a limit may lead to overfitting noise as small bandwidth peaks.
        We recommend a lower bound of approximately 2x the frequency resolution.

Model Parameters

Once the power spectrum model has been calculated, the model fit parameters are stored as object attributes that can be accessed after fitting.

Access model fit parameters from specparam object, after fitting:

# Aperiodic parameters
print('Aperiodic parameters: \n', fm.results.params.aperiodic.params, '\n')

# Peak parameters
print('Peak parameters: \n', fm.results.params.periodic.params, '\n')

# Check how many peaks were fit
print('Number of fit peaks: \n', fm.results.n_peaks)

Aperiodic parameters:
 [-21.61849957   0.7160227 ]

Peak parameters:
 [[ 9.36187382  1.04444622  1.58733674]
 [11.1723471   0.23008765  2.87693709]
 [18.24843214  0.33140755  2.84640768]]

Number of fit peaks:
 3

Selecting Parameters

You can also select parameters using the get_params() method, which can be used to specify which parameters you want to extract.

# Extract parameters, indicating sub-selections of parameters
exp = fm.get_params('aperiodic', 'exponent')
cfs = fm.get_params('periodic', 'CF')

For a full description of how you can access data with get_params(), check the method’s documentation.

As a reminder, you can access the documentation for a function using ‘?’ in a Jupyter notebook (ex: fm.get_params?), or more generally with the help function in general Python (ex: help(fm.get_params)).

Model Metrics

In addition to model fit parameters, the fitting procedure computes and stores various metrics that can be used to evaluate model fit quality.

# Goodness of fit metrics
print('Goodness of fit:')
print(' Error - ', fm.results.metrics.results['error_mae'])
print(' R^2   - ', fm.results.metrics.results['gof_rsquared'], '\n')

Goodness of fit:
 Error -  0.035608193181702164
 R^2   -  0.9828921187440269

You can also access metrics with the get_metrics() method.

# Extract a model metric with `get_metrics`
err = fm.get_metrics('error')

Extracting model fit parameters and model metrics can also be combined to evaluate model properties, for example using the following template:

# Print out a custom parameter report
template = ("With an error level of {error:1.2f}, an exponent "
            "of {exponent:1.2f} and peaks of {cfs:s} Hz were fit.")
print(template.format(error=err, exponent=exp,
                      cfs=' & '.join(map(str, [round(cf, 2) for cf in cfs]))))

With an error level of 0.04, an exponent of 0.72 and peaks of 9.36 & 11.17 & 18.25 Hz were fit.

# Compare the 'peak_params_' to the underlying gaussian parameters
print('  Peak Parameters \t Gaussian Parameters')
for peak, gauss in zip(fm.results.params.periodic.get_params('converted'),
                       fm.results.params.periodic.get_params('fit')):
    print('{:5.2f} {:5.2f} {:5.2f} \t {:5.2f} {:5.2f} {:5.2f}'.format(*peak, *gauss))

  Peak Parameters        Gaussian Parameters
 9.36  1.04  1.59         9.36  0.98  0.79
11.17  0.23  2.88        11.17  0.17  1.44
18.25  0.33  2.85        18.25  0.33  1.42

FitResults

There is also a convenience method to return all model fit results: get_results().

This method returns all the model fit parameters, including the underlying Gaussian parameters, collected together into a FitResults object.

The FitResults object, which in Python terms is a named tuple, is a standard data object used to organize and collect parameter data.

# Grab each model fit result with `get_results` to gather all results together
#   Note that this returns a FitResults object
fres = fm.results.get_results()

# You can also unpack all fit parameters when using `get_results`
ap_fit, ap_conv, peak_fit, peak_conv, metrics = fm.results.get_results()

# Print out the FitResults
print(fres, '\n')

# from specparamResults, you can access the different results
print('Aperiodic Parameters: \n', fres.aperiodic_fit)

# Check the R^2 and error of the model fit
print('R-squared: \n {:5.4f}'.format(fres.metrics['gof_rsquared']))
print('Fit error: \n {:5.4f}'.format(fres.metrics['error_mae']))

FitResults(aperiodic_fit=array([-21.61849957,   0.7160227 ]), aperiodic_converted=array([nan, nan]), peak_fit=array([[ 9.36187382,  0.97903265,  0.79366837],
       [11.1723471 ,  0.16897043,  1.43846854],
       [18.24843214,  0.3341275 ,  1.42320384]]), peak_converted=array([[ 9.36187382,  1.04444622,  1.58733674],
       [11.1723471 ,  0.23008765,  2.87693709],
       [18.24843214,  0.33140755,  2.84640768]]), metrics={'error_mae': 0.035608193181702164, 'gof_rsquared': 0.9828921187440269})

Aperiodic Parameters:
 [-21.61849957   0.7160227 ]
R-squared:
 0.9829
Fit error:
 0.0356

Conclusion

In this tutorial, we have explored the basics of the SpectralModel object, fitting power spectrum models, and extracting parameters.

In the next tutorial, we will explore how this algorithm actually works to fit the model.