diff --git a/doc/api.rst b/doc/api.rst index 33f16f1d..b7eba362 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -67,16 +67,18 @@ Functions to manipulate, examine, and analyze model objects. average_reconstructions fit_models_3d +Component Objects +----------------- -Data Objects ------------- +The model object combines multiple sub-objects that define and store different +elements of the model definition, data, and results. -Objects to manage frequency bands, model information, and simulation parameters. +Here, the main sub-objects are listed, some of which can also be used independently. -Bands Object -~~~~~~~~~~~~ +Bands +~~~~~ -An object to handle frequency band definitions. +An object for defining frequency band definitions. .. currentmodule:: specparam @@ -85,6 +87,98 @@ An object to handle frequency band definitions. Bands +Algorithm +~~~~~~~~~ + +An object for defining fit algorithms. + +.. currentmodule:: specparam.algorithms.algorithm + +.. autosummary:: + :toctree: generated/ + + Algorithm + +Modes +~~~~~ + +An object for defining fit modes. + +.. currentmodule:: specparam.modes.mode + +.. autosummary:: + :toctree: generated/ + + Mode + +Associated objects allow for defining mode parameters. + +.. currentmodule:: specparam.modes.params + +.. autosummary:: + :toctree: generated/ + + ParamDefinition + +Utility to check for available fit modes. + +.. currentmodule:: specparam.modes.definitions + +.. autosummary:: + :toctree: generated/ + + check_modes + +Metrics +~~~~~~~ + +An object for defining metrics. + +.. currentmodule:: specparam.metrics.metric + +.. autosummary:: + :toctree: generated/ + + Metric + +Utility to check for available metrics. + +.. currentmodule:: specparam.metrics.definitions + +.. autosummary:: + :toctree: generated/ + + check_metrics + +Data +~~~~ + +An object for managing data to be modeled. + +.. currentmodule:: specparam.data.data + +.. autosummary:: + :toctree: generated/ + + Data + +Results +~~~~~~~ + +An object for managing model results. + +.. currentmodule:: specparam.results.results + +.. autosummary:: + :toctree: generated/ + + Results + +Data Objects +------------ + +Objects to manage model information, and simulation parameters. + Model Information ~~~~~~~~~~~~~~~~~ diff --git a/doc/conf.py b/doc/conf.py index 46491058..c4d2a7dc 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -4,6 +4,7 @@ # # For a full list of documentation options, see: # https://www.sphinx-doc.org/en/master/usage/configuration.html +# # ---------------------------------------------------------------------------- @@ -80,6 +81,7 @@ # Settings for sphinx_copybutton copybutton_prompt_text = "$ " + # -- Options for HTML output ------------------------------------------------- # The theme to use for HTML and HTML Help pages. @@ -126,6 +128,7 @@ def setup(app): app.add_css_file("my-styles.css") + # -- Extension configuration ------------------------------------------------- # Configurations for sphinx gallery @@ -136,6 +139,7 @@ def setup(app): '../examples/manage', '../examples/models', '../examples/plots', + '../examples/customize', '../examples/sims', '../examples/analyses', '../motivations/concepts', diff --git a/doc/faq.rst b/doc/faq.rst index 052e9ccc..d8ec7db5 100644 --- a/doc/faq.rst +++ b/doc/faq.rst @@ -3,9 +3,8 @@ Frequently Asked Questions The following is a collection of frequently asked questions and answers about spectral parameterization. -These answers focus on the ideas and concepts relating to parameterizing neural power spectra. - -For code related questions, check out the API listing, tutorials, and examples. +The questions and answers here focus on the key ideas and concept relating to the approach. +For code related questions and examples, check out the API listing and tutorials. Table of Contents ----------------- @@ -16,23 +15,46 @@ Table of Contents What is spectral parameterization? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Spectral parameterization means fitting a model to describe power spectra. -A particular algorithm and implementation for doing this is available in `specparam`, -an open-source Python module for parameterizing neural power spectra. +Spectral parameterization means fitting a model to describe a power spectrum. +The `specparam` module implements a framework for fitting models to power spectra, +available in an open-source Python module. -Spectral parameterization uses a model-driven approach that assumes that neurophysiological time -series are comprised of two separable components, reflecting periodic (or oscillatory) and -aperiodic activity. This approach relies on the assumption that these two components are indeed separable -components of the underlying data, though it is agnostic to their physiological origin and -putative functional roles. +How does spectral parameterization work? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Spectral parameterization operates on frequency representations of neurophysiological times -series (power spectra). At it's core, the module contains an algorithm to measure these two -components - the periodic and aperiodic components - in power spectra. The final model -of the neural power spectrum consists of quantifications of each of the two components, as well as -a combined model fit of the whole power spectrum. +Spectral parameterization operates on frequency representations (power spectra) of digital +signals (time series). Broadly speaking, the module contains model definitions for how to +characterize each component (periodic and aperiodic) and algorithms to fit these definitions +to power spectra. The resulting model of the power spectrum consists of quantifications of +each of the two components, as well as a combined model fit of the whole power spectrum. +For more information on the approach, including descriptions of the model definitions +and algorithms, see the Tutorials. + +What was spectral parameterization originally designed for? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The original use case of spectral parameterization was for analyzing +neuro-electrophysiological data. Due to this, across the documentation, +the vast majority of the discussion of the use of spectral parameterization +refers to the use case of neuro-electrophysiological data. + +What data can this be applied to? Can it be used on non-neural data? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Spectral parameterization can, in theory, be applied to any power spectra, +including non-neural power spectra. Following developments and generalizations of the module, +starting with the specparam 2.0 release version, the module should now be more easily +customizable and applicable to other data types. In practice, any application of spectral +parameterization needs to evaluate if the model fitting procedures are well posed +for the specific use case (see Tutorials for more details). + +What is needed to parameterize a power spectrum? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The full mathematical description of the model is described in the tutorials. +In order for spectral parameterization to work for a given application, the chosen +model definition has to be appropriate for the data under study, and the chosen fitting +algorithm has to be appropriate for fitting this model to the data. See discussions throughout +the tutorials on choosing model forms and fit algorithms. What is meant by 'aperiodic' activity? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -58,12 +80,26 @@ and detectable by the fitting algorithm if they exhibit as band-limited power ov the aperiodic component in the power spectrum. This 'peak' of power over the aperiodic is taken as evidence of frequency specific power, distinct from the aperiodic component. -Note that this periodic activity need not be continuous, as oscillatory activity often -exhibits as 'bursts' in the time series, nor sinusoidal, as rhythmic neural activity is -often non-sinusoidal. +Note that, in the time domain, this periodic activity need not strictly be continuous nor +completely sinusoidal, as variability in oscillatory activity is still typically reflected +as peaks of power in the spectral domain. -Why should I parameterize power spectra? -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +For example, in neural data, oscillatory activity often exhibits as 'bursts' in the time series, +which are often at least somewhat non-sinusoidal, which can still be seen as peaks in the power +spectrum (though see notes on interpreting peaks of power). + +What are the assumptions of spectral parameterization? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Spectral parameterization uses a model-driven approach that assumes that the data under study +are comprised of two separately measure-able components, reflecting periodic (or oscillatory) +and aperiodic activity. This approach therefore relies on the assumption that these two components +are indeed separable components of the power spectrum as observed by the model. The model is +broadly agnostic to the relationship(s) between the components, their origin(s) in the generative +process that created the data, and their putative functional roles. + +Why should I parameterize neural power spectra? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Though research often focuses on periodic (rhythmic or oscillatory) components, neurophysiological recordings of electrophysiological neural activity contain not only periodic, but also aperiodic @@ -90,8 +126,8 @@ For more discussion and examples of the conceptual and methodological factors th motivate parameterizing neural power spectra, see the `motivations page `_ -Why is it important to measure aperiodic activity? -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Why is it important to measure aperiodic neural activity? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Aperiodic activity has long known to be present in neural data, but has been less of a research focus (as compared to periodic activity). Recent work has demonstrated @@ -207,20 +243,22 @@ is a Python tool for analyzing neural oscillations and their waveform shape prop Why is parameterizing neural power spectra different from other methods? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -There are many existing methods for analyzing periodic activity, and also other methods for -analyzing aperiodic activity. Most existing methods are designed to measure one or the other -signal component. Few methods attempt to explicitly separate and quantify both the periodic -and aperiodic components of the signal. This combined approach is a key factor that we -consider to be important for getting the measurements to work well. By jointly measuring -both components, the method is more capable of quantifying which aspects of the data -are changing and in what ways. +Within neuroscience, there are many existing methods for analyzing periodic activity, +and also other methods for analyzing aperiodic activity. Historically, fewer methods have +been developed that attempt to explicitly separate and quantify both the periodic and +aperiodic components of the signal. As such, at the time of the original development of specparam, +most existing methods were designed to measure one pre-specified signal component +(periodic or aperiodic). Using the combined approach of spectral parameterization (considering +and measuring both components together) is therefore a key factor that is different from +many other approaches. The goal is that by jointly measuring both components, the method should be +more capable of quantifying which aspects of the data are changing and in what ways. -More in depth analyses of the properties of the fitting algorithm, and systematic comparisons -with other methods (through simulations) are are also ongoing, to clarify when and how -this approach compares to different methods. +Since the original development of the method, there has been significant development of other +methods (partially summarized here: [7_]) and work comparing between different methods. +Check the current literature for more discussion of these topics. -What data modalities can be used for parameterizing neural power spectra? -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +What neuroscience modalities can be analyzed with spectral parameterization? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The power spectrum model can theoretically be applied to power spectra derived from any electrophysiological or magnetophysiological signal of neural origin. In practice, this @@ -369,15 +407,18 @@ one or the other might be more appropriate. Is there a time resolved version? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Since it operates on frequency representations (power spectra) the power spectrum model is not, -by construction, a time resolved method. +Yes, as of the specparam 2.0 version, the module includes functionality to fit models across +time (across spectrograms). -However, similar to other frequency estimation approaches that are used in a time-resolved manner, -it can, in theory, be applied in a sliding window fashion. This approach could be used to estimate -spectral features across time, somewhat analogous to a spectrogram. +Note that as it operates in the frequency domain, in order to be able to analyze data over time, +the model is applied to individual windows over data, whereby each window reflects a time-segment +of data. By fitting models across a series of windows, spectral parameterization results +can be examined across time (across windows). In doing so, it is therefore important to +consider the spectral estimation, in terms of key details such as window size, method of +estimation, window overlap, etc, in order to make sure the models are well fit and to +appropriately interpret the results. -This functionality is not currently available or described in the current module, but is a focus -off current work. We hope to add information, guidelines, and tooling to do this once this soon. +See more information on this functionality in the tutorials. What is the 'FOOOF' name? ~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -389,6 +430,13 @@ we have moved away from using these terms in the module and algorithm, now prefe as 'periodic' and 'aperiodic' activity, the module has been renamed to the more general name of 'spectral parameterization'. +In addition, the generalizability of the approach has changed with the different versions of the +code. The `fooof` versions of the module can be though of as implementing a particular model +definition and specific algorithm for parameterizing neural power spectra. By comparison, the +`specparam` versions of the module are designed to enable flexible model definition and application +to power spectra across different contexts (while still also including the original +specific functionality from fooof). + How do I cite this method? ~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -424,3 +472,8 @@ References - [6_] He (2014). Scale-free Brain Activity: Past, Present and Future. DOI: 10.1016/j.tics.2014.04.003 .. _6 : https://doi.org/10.1016/j.tics.2014.04.003 + +- [7_] Donoghue & Watrous (2023). How Can We Differentiate Narrow-Band Oscillations from Aperiodic Activity? + DOI: 10.1007/978-3-031-20910-9_22 + +.. _7 : https://doi.org/10.1007/978-3-031-20910-9_22 diff --git a/examples/customize/README.txt b/examples/customize/README.txt new file mode 100644 index 00000000..45690cf5 --- /dev/null +++ b/examples/customize/README.txt @@ -0,0 +1,4 @@ +Customize +--------- + +Examples of customizing the fitting processes. \ No newline at end of file diff --git a/examples/customize/plot_custom_metrics.py b/examples/customize/plot_custom_metrics.py new file mode 100644 index 00000000..1535d29e --- /dev/null +++ b/examples/customize/plot_custom_metrics.py @@ -0,0 +1,340 @@ +""" +Custom Metrics +============== + +This example covers defining and using custom metrics. +""" + +# Import spectral model object +from specparam import SpectralModel + +# Import Metric object for defining custom metrics +from specparam.metrics.metric import Metric + +# Import function to simulate power spectra +from specparam.sim import sim_power_spectrum + +################################################################################################### +# Defining Custom Metrics +# ----------------------- +# +# As covered in the tutorials, the specparam module has a set of predefined metrics, wherein +# `metrics` refer to measures that are computed post model fitting to evaluate properties +# of the model fit, typically as it relates to the original data. +# +# In this tutorial, we will explore how you can also define your own custom metrics. +# +# To do so, we will start by simulating an example power spectrum to use for this example. +# + +################################################################################################### + +# Define simulation parameters +ap_params = [50, 2] +gauss_params = [10, 0.5, 2, 20, 0.3, 4] +nlv = 0.05 + +# Simulate an example power spectrum +freqs, powers = sim_power_spectrum([3, 50], {'fixed' : ap_params}, {'gaussian' : gauss_params}, + nlv, freq_res=0.25) + +################################################################################################### +# Defining Custom Measure Function +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# To define a custom measure, we first need to define a function that computes our metric. +# +# By convention, a custom measure functions must have at least two input arguments, with the +# first two arguments being the original data (power spectrum) and the model (modeled spectrum). +# +# Within the function, the measure of interest should be defined, such that it returns +# the result of the metric, which should be a float. +# +# For our first example, we will define a simple error metric that computes the total +# error of the model (sum of the absolute deviations). +# + +################################################################################################### + +# Import any functionality needed for the metric +import numpy as np + +# Define a function that computes a custom metric +def compute_total_error(power_spectrum, modeled_spectrum): + """Compute the total (summed) error between the data and the model.""" + + total_err = np.sum(np.abs(power_spectrum - modeled_spectrum)) + + return total_err + +################################################################################################### +# The `Metric` Class +# ~~~~~~~~~~~~~~~~~~ +# +# In order to use the custom metric, more information is needed. To collect this additional +# information the :class:`~specparam.metrics.metric` object is used to define a metric. +# +# The Metric object requires the following information: +# +# - `category`: a description of what kind of metric it is +# - `measure`: a label for the specific measure that is defined +# - `description`: a description of the custom metric +# - `func`: the callable that compute the metric +# + +################################################################################################### + +# Define Metric for the total error +total_error_metric = Metric( + category='error', + measure='total_error', + description='Total absolute error.', + func=compute_total_error, +) + +################################################################################################### + +# Initialize a spectral model, passing in our custom metric definition +fm = SpectralModel(min_peak_height=0.25, metrics=[total_error_metric]) + +# Fit the model and print a report +fm.report(freqs, powers) + +################################################################################################### +# +# Note that in the above report, the metrics section now includes the result of our custom metric! +# + +################################################################################################### +# Defining Metrics with Dictionaries +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# In the above, we directly used the `Metric` object to define our custom metric. +# +# If you prefer, you can also collect the relevant information needed to define a metric into +# a dictionary, and pass this into the model object instead. +# + +################################################################################################### + +# Define the information for a custom metric, in a dictionary +custom_metric_dict = { + 'category' : 'error', + 'measure' : 'total_error', + 'description' : 'Total absolute error.', + 'func' : compute_total_error, +} + +################################################################################################### + +# Initialize a model object, passing in the custom metric, defined as a dictionary +fm = SpectralModel(min_peak_height=0.25, metrics=[custom_metric_dict]) + +################################################################################################### +# +# When using custom metrics, you can also access the results using the +# :meth:`~specparam.SpectralModel.get_metrics` by using the measure name, +# the same as when accessing default / built in metrics. +# + +################################################################################################### + +# Fit the model to the data +fm.fit(freqs, powers) + +# Access the custom metric result with get_metrics +fm.get_metrics('total_error') + +################################################################################################### +# +# Above, we initialized our model object by specifying to use only our new custom metric. +# +# Note that you can also initialize the model object with a list of multiple metrics, +# including a mixture of pre-defined and/or custom defined metrics. +# + +################################################################################################### + +# Initialize a spectral model, passing in multiple metrics (both internal and custom) +fm = SpectralModel(min_peak_height=0.25, metrics=[total_error_metric, 'gof_rsquared']) + +# Fit the model and print a report +fm.report(freqs, powers) + +################################################################################################### +# +# In the above report, we now see multiple metrics have been applied to the model, including +# our new / custom metric as well as the built-in metrics which we specified. +# + +################################################################################################### +# Custom Metrics with Additional Arguments +# ---------------------------------------- +# +# In some cases, you may want to define a metric that requires additional information than just +# the data and model to compute the measure of interest (e.g., the custom metric function has +# more than two arguments). +# +# For an example of this, we will define an example custom metric that computes the error +# of a specific frequency range, therefore requiring information about the frequency definition +# of the data. +# +# We start, as above, by defining a function that computes our metric, starting with the same +# two arguments (data and model) and then adding additional arguments as needed. +# + +################################################################################################### + +from specparam.utils.spectral import trim_spectrum + +def compute_lowfreq_error(power_spectrum, modeled_spectrum, freqs): + """Compute the mean absolute error in the low frequency range.""" + + low_freq_range = [1, 8] + _, power_spectrum_low = trim_spectrum(freqs, power_spectrum, low_freq_range) + _, modeled_spectrum_low = trim_spectrum(freqs, modeled_spectrum, low_freq_range) + + low_err = np.abs(power_spectrum_low - modeled_spectrum_low).mean() + + return low_err + +################################################################################################### +# +# In the above error function, we need access to the frequency definition, as well as the +# data and model. Now we need to make sure that when this function is called to compute +# the metric, this additional information is made available to the function. +# +# To provide access to additional attributes, we need to use the optional `kwargs` argument +# when defining our Metric to define how to access the additional information. +# +# To use kwargs, define a dictionary where each key is the string name of the additional input +# to the measure function, and each value is a lambda function that accepts as input the model +# data & results objects, and then uses these to access the information that is needed. +# +# For our low frequency error measure, this looks like: +# +# `kwargs={'freqs' : lambda data, results: data.freqs}` +# +# Internally, if `kwargs` is defined, the lambda function is called for each entry, +# passing in the Model.data and Model.results objects, which then accesses the specified +# attributes based on the implementation of the lambda function. +# +# Note that this means all additional inputs to the function need to be information that +# can be accessed and/or computed based on what is available in the data and results +# objects that are part of the model object +# + +################################################################################################### + +# Define Metric for the low frequency error, defining `kwargs` +lowfreq_error = Metric( + category='error', + measure='low_freq_mae', + description='Mean absolute error of the low frequency range.', + func=compute_lowfreq_error, + kwargs={'freqs' : lambda data, results: data.freqs}, +) + +################################################################################################### +# +# Now our custom metric is defined, and we can use it with a model object the same as before! +# + +################################################################################################### + +# Initialize a spectral model, passing in custom metric with additional arguments +fm = SpectralModel(metrics=[lowfreq_error]) + +# Fit the model and print a report +fm.report(freqs, powers) + +################################################################################################### +# +# In the above, our new custom metric was now computed for our model fit! +# + +################################################################################################### +# A Final Example +# --------------- +# +# For one last example, lets make a more complex metric, which requires multiple additional +# pieces of information that need to be accessed from the model object. +# +# In this example, we will define and use a custom metric that defines an error metric +# that is proportional to the model degrees of freedom. +# + +################################################################################################### + +# Define a function to compute our custom error metric +def custom_measure(power_spectrum, modeled_spectrum, freqs, n_params): + """Compute a custom error metric of error proportional to model degrees of freedom.""" + + # Compute degrees of freedom (# data points - # parameters) + df_error = len(freqs) - n_params + + # Compute the total error of the model fit + err_total = compute_total_error(power_spectrum, modeled_spectrum) + + # Compute the error proportional to model degrees of freedom + err_per_df = err_total / df_error + + return err_per_df + +################################################################################################### +# +# Now that we have defined the function, we define a Metric object, as before. +# +# Note that in defining the 'category' for our custom metric, we need not use the existing +# categories from the built in metrics, and can instead define our own custom category. +# + +################################################################################################### + +# Define Metric for the low frequency error +custom_measure = Metric( + category='custom', + measure='err-by-df', + description='Error proportionate to the degrees of freedom of the model.', + func=custom_measure, + kwargs={'freqs' : lambda data, results: data.freqs, + 'n_params' : lambda data, results : results.n_params}, +) + +################################################################################################### + +# Initialize a spectral model, passing in our new custom measure +fm = SpectralModel(metrics=[custom_measure]) + +# Fit the model and print a report +fm.report(freqs, powers) + +################################################################################################### +# +# We can again use `get_metrics` to access the metric results - note that in this case +# we need to match the category name that we used in defining our metric. +# + +################################################################################################### + +# Access the custom metric result +fm.get_metrics('custom_err-by-df') + +################################################################################################### +# +# That covers how to define custom metrics! + +################################################################################################### +# Adding New Metrics to the Module +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# As a final note, if you look into the set of 'built-in' metrics that are available within +# the module, you will see that these are defined in the exact way as done here - the only +# difference is that they are defined within the module and therefore can be accessed via +# their name, as a shortcut, rather than the user having to pass in their own full definitions. +# +# This also means that if you have a custom metric that you think would be of interest to +# other specparam users, once the Metric object is defined it is quite easy to add this +# to the module as a new default option. If you would be interested in suggesting a metric +# be added to the module, feel free to open an issue and/or pull request. +# diff --git a/examples/customize/plot_custom_modes.py b/examples/customize/plot_custom_modes.py new file mode 100644 index 00000000..5d1204cd --- /dev/null +++ b/examples/customize/plot_custom_modes.py @@ -0,0 +1,339 @@ +""" +Custom Modes +============ + +This example covers defining and using custom fit modes. +""" + +# sphinx_gallery_thumbnail_number = 3 + +# Import functionality used to define custom fit modes +from collections import OrderedDict +import numpy as np + +# Import the model object +from specparam import SpectralModel + +# Import function to simulate a power spectrum +from specparam.sim import sim_power_spectrum + +# Import objects used to define a custom mode +from specparam.modes.mode import Mode +from specparam.modes.params import ParamDefinition + +################################################################################################### +# Defining Custom Fit Modes +# ------------------------- +# +# As covered in the tutorials, the specparam module has a set of predefined fit modes, wherein +# `modes` refer to the fit function and associated metadata & description that describes how +# each component is fit. +# +# In this tutorial, we will explore how you can also define your own custom modes. +# +# To do so, we will start by simulating an example power spectrum to use for this example. +# + +################################################################################################### + +# Define simulation parameters +ap_params = [0, 1] +gauss_params = [10, 0.5, 2, 20, 0.3, 2] +nlv = 0.0025 + +# Simulate an example power spectrum +freqs, powers = sim_power_spectrum(\ + [1, 40], {'fixed' : ap_params}, {'gaussian' : gauss_params}, nlv, freq_res=0.25) + +################################################################################################### +# Example: Custom Aperiodic Mode +# ------------------------------ +# +# To start, we will define a custom aperiodic fit mode. Specifically, we will define a +# custom aperiodic fit mode that fits only a single exponent. +# +# In theory, such an aperiodic fit mode could be used if one knew that the offset for all +# power spectra of interest had an offset of 0 (for example, if they were all normalized as such), +# but in practice this is not likely to be a useful fit mode, and is used here merely as +# a simplified example of defining a custom fit mode. That is, while this is a functioning +# custom fit mode definition for the sake of an example, this is not necessarily a useful +# (or recommended) fit function. +# +# First, we need to define a fit function that can be applied to fit our component. +# +# To do so, we define a function that defines the fit function. To be able to be used in a +# model object, this function needs to follow a couple key principles: +# +# - the first input should be for x-axis values of the data to be modeled (frequencies) +# - regardless of how many parameters the fit function has, these should be organized as *args +# in the function definition +# - the body of function should unpack the params, fit the function definition to the input +# data, and then return the model result (modeled component) +# +# By following these conventions, the model object is able to use this function to fit the +# specified component, passing in data and parameters appropriately, regardless of the +# set of parameters and function definition. +# + +################################################################################################### + +# Define a custom aperiodic fit function +def expo_only_function(xs, *params): + """Exponent only aperiodic fit function + + Parameters + ---------- + xs : 1d array + X-axis values. + *params : float + Parameters that define the component fit. + + Returns + ------- + ys : 1d array + Output values of the fit. + """ + + exp = params + ys = np.log10(xs**exp) + + return ys + +################################################################################################### +# +# Now that we have the fit function defined, we need to collect some additional information +# to be able to use our fit mode in the model object, starting with the parameter definition. +# +# In order for the model object to have a description of the parameters that define the fit +# mode, we use the :class:`~specparam.modes.params.ParamDefinition` object. +# +# To use this object, we use an OrderedDict to define a parameter description, where each key +# is a parameter name (this being the name that will be used to access the parameter from the +# model object), and a description of the parameter. Note that by using an OrderedDict, we +# ensure that the order of the parameters is consistent. Make sure the order of the parameters +# in the definition matches the order of the parameters in the fit function. +# + +################################################################################################### + +# Define the parameter definition for the expo only fit mode +expo_only_params = ParamDefinition(OrderedDict({ + 'expo' : 'Exponent of the aperiodic component.', +})) + +################################################################################################### +# +# Now we have the main elements of our new fit mode, we can use the +# :class:`~specparam.modes.mode.Mode` object to define the full mode. To do so, we initialize +# a Mode object passing in metadata about the fit mode, as well our our fit function and +# parameter definition. +# +# Note that there is some meta-data that needs to be defined in the Mode definition, including: +# +# - name, component, description: describes the fit mode +# - jacobian: a function that computes the Jacobian for the fit function, +# which in some cases can speed up fitting. +# - ndim: the dimensionality of the output parameters, which should typically be 1 for +# aperiodic modes (returns a 1d array of parameters per model fit) and 2d for peak parameters +# (returns a 2d array of parameters across potentially multiple detected peaks) +# - freqs_space, powers_space : the expected spacing of the data +# + +################################################################################################### + +# Define the custom exponent only fit mode +expo_only_mode = Mode( + name='custom_expo_only', + component='aperiodic', + description='Fit an exponent only.', + func=expo_only_function, + jacobian=None, + params=expo_only_params, + ndim=1, + freq_space='linear', + powers_space='log10', +) + +################################################################################################### +# +# Our custom fit mode if now defined! +# +# The use this fit mode, we initialize a model object, passing the fit mode in as the specified +# component mode. +# + +################################################################################################### + +# Initialize model object, passing in custom aperiodic mode definition +fm = SpectralModel(aperiodic_mode=expo_only_mode) + +################################################################################################### +# +# Now we can use the model object to fit some data. +# + +################################################################################################### + +# Fit model and report results +fm.report(freqs, powers) + +################################################################################################### +# +# In the above report we can see that under the aperiodic mode section, the results reflect +# our custom fit mode! +# +# Note that the parameter name is taken from what we described in the parameter definition. This +# is also the name you use to access the parameter results in `get_params`. +# + +################################################################################################### + +# Get the aperiodic parameters, using the custom fit mode parameter name +fm.get_params('aperiodic', 'expo') + +################################################################################################### +# Example Custom Periodic Fit Mode +# -------------------------------- +# +# Defining a custom fit mode can also be done in the same way for periodic modes. +# +# In this example, we will define and apply a custom periodic fit mode that uses a +# rectangular peak fit function. Note that as we will see, this is not really a good +# peak option for neural data (though it may be useful for other data), but serves here +# as an example of +# +# First, we start by defining a fit function, as before. This should follow the same format +# as introduced previously for the aperiodic fit mode function, with the added consideration +# that peak functions should be flexible for potentially having multiple possible peaks, +# meaning that the fit function needs to be set up in a way that multiple sets of parameters +# for multiple peaks can be passed in together, and the results summed together (e.g., if the +# fit function takes `p` parameters, the function should accept multiples of `p` number of inputs +# and loops across these, summing the resultant peaks). +# + +################################################################################################### + +# Define a custom peak fit function +def rectangular_function(xs, *params): + """Custom periodic fit function - rectangular fit. + + Parameters + ---------- + xs : 1d array + X-axis values. + *params : float + Parameters that define the component fit. + + Returns + ------- + ys : 1d array + Output values of the fit. + """ + + ys = np.zeros_like(xs) + + for ctr, hgt, wid in zip(*[iter(params)] * 3): + ys[np.abs(xs - ctr) <= wid] += 1 * hgt + + return ys + +################################################################################################### +# +# Next up, as before, we need to define a parameter definition. Here, we will use the same labels +# as the standard (Gaussian) peak mode, redefined for the rectangle case. +# + +################################################################################################### + +rectangular_params = ParamDefinition(OrderedDict({ + 'cf' : 'Center frequency of the rectangle.', + 'pw' : 'Power of the rectangle, over and above the aperiodic component.', + 'bw' : 'Width of the rectangle.' +})) + +################################################################################################### +# +# Finally, as before, we collect together all the required information to define our fit mode. +# + +################################################################################################### + +# Define the custom rectangular fit mode +rectangular_mode = Mode( + name='rectangular', + component='periodic', + description='Custom rectangular (boxcar) peak fit function.', + func=rectangular_function, + jacobian=None, + params=rectangular_params, + ndim=2, + freq_space='linear', + powers_space='log10', +) + +################################################################################################### +# +# Now we can initialize a model object with our custom periodic fit mode, and use it to fit +# some data! +# + +################################################################################################### + +# Initialize model object, passing in custom periodic mode definition +fm = SpectralModel(periodic_mode=rectangular_mode, max_n_peaks=2) + +################################################################################################### + +# Fit model and report results +fm.report(freqs, powers) + +################################################################################################### +# +# In the above, we can see the results of using the custom periodic mode peak fit function. +# + +################################################################################################### +# Relationship Between Fit Modes & Fit Algorithms +# ----------------------------------------------- +# +# In these examples, we defined simple fit modes wherein the new functions are similar enough +# to the existing cases that they can be dropped in and still work with the default fit +# algorithm. There are some new fit modes that might not work well with the existing / default +# algorithm - for example, the default algorithm makes some assumptions about the peak fit +# function. This means that major changes to the fit modes may need to be defined +# in tandem with updates to the fitting algorithm to make it all work together. Relatedly, +# you might notice from the above mode definition that there are additional details +# that can be customized, such as the expected spacing of the data, that would also require +# coordination with making sure the fit algorithm is consistent with the defined fit modes. +# + +################################################################################################### +# Combining Custom Modes +# ---------------------- +# +# As a final example, note that you can also combine custom periodic and aperiodic fit modes +# together. +# + +################################################################################################### + +# Initialize model object, passing in custom aperiodic & periodic mode definitions +fm = SpectralModel(aperiodic_mode=expo_only_mode, periodic_mode=rectangular_mode, max_n_peaks=2) + +# Fit model and report results +fm.report(freqs, powers) + +################################################################################################### +# Adding New Modes to the Module +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# As a final note, if you look into the set of 'built-in' modes that are available within +# the module, you will see that these are defined in the exact way as done here - the only +# difference is that they are defined within the module and therefore can be accessed via +# their name, as a shortcut, rather than the user having to pass in their own full definitions. +# +# This also means that if you have a custom mode that you think would be of interest to +# other specparam users, once the Mode object is defined it is quite easy to add this +# to the module as a new default option. If you would be interested in suggesting a mode +# be added to the module, feel free to open an issue and/or pull request. +# diff --git a/examples/manage/plot_manipulating_models.py b/examples/manage/plot_manipulating_models.py index eb9214fb..02eaa692 100644 --- a/examples/manage/plot_manipulating_models.py +++ b/examples/manage/plot_manipulating_models.py @@ -186,12 +186,12 @@ ################################################################################################### # Check information on null models (dropped models) -print('Number of Null models : \t', fg.results.n_null) -print('Indices of Null models : \t', fg.results.null_inds) +print('Number of Null models :\t', fg.results.n_null) +print('Indices of Null models :\t', fg.results.null_inds) # Despite the dropped model, the total number of models in the object is the same # This means that the indices are still the same as before dropping models -print('Number of model fits: ', len(fg.results)) +print('Number of model fits :\t', len(fg.results)) ################################################################################################### diff --git a/tutorials/plot_05-AperiodicFitting.py b/tutorials/plot_05-AperiodicFitting.py index 3f67fbdd..8cd663a7 100644 --- a/tutorials/plot_05-AperiodicFitting.py +++ b/tutorials/plot_05-AperiodicFitting.py @@ -163,7 +163,7 @@ ################################################################################################### # A note on interpreting the 'knee' parameter -# ------------------------------------------- +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # The aperiodic fit has the form: # @@ -214,6 +214,35 @@ # 'knee' mode to be able to appropriately characterize the data. # +################################################################################################### +# Aperiodic Fit Mode: doublexp +# ---------------------------- +# +# Returning to our exploration of the available fit modes for the aperiodic component, +# another avaible fit mode is the 'double exponent' or 'doublexp'. +# +# In the above 'knee' mode, you might have noticed that even though the knee is described as +# a change in the value of the aperiodic exponent, implying there are two different +# exponent values, we still only fit one exponent value. In the above case, the exponent +# above the knee is fit, whereas the exponent below the knee is assumed to be zero. +# +# By comparison, the double exponent model fits two different exponent values, above and +# below the knee. +# + +################################################################################################### + +# Create and fit a power spectrum model in doublexp aperiodic fit mode +fm = SpectralModel(peak_width_limits=[2, 8], aperiodic_mode='doublexp') +fm.report(freqs, spectrum, [2, 70], plt_log=True) + +################################################################################################### +# +# In the above example model fit, we can see that the aperiodic mode now includes 4 fit parameters, +# including two different exponent values (exponent1, reflecting below the knee, and exponent2, +# reflecting above the knee), as well as the offset and knee value. +# + ################################################################################################### # Choosing an Aperiodic Fit Mode # ------------------------------ @@ -233,7 +262,7 @@ # # - This is likely across smaller frequency ranges, such as 3-30. # - Do not perform no-knee fits across a range in which this does not hold. -# - If there is a clear knee, then use knee fits. +# - If there is a clear knee, then use a fit mode that includes a knee. # # - This is likely across larger fitting ranges such as 1-150 Hz. # - Be wary of ambiguous ranges, where there may or may not be a knee. diff --git a/tutorials/plot_07-ModelObject.py b/tutorials/plot_07-ModelObject.py index ed1dcac8..002a7598 100644 --- a/tutorials/plot_07-ModelObject.py +++ b/tutorials/plot_07-ModelObject.py @@ -1,5 +1,5 @@ """ -07: Exploring the model object +07: Exploring the Model Object ============================== Further exploring the SpectralModel object, including algorithm settings and available methods. diff --git a/tutorials/plot_08-GroupFits.py b/tutorials/plot_08-GroupFits.py index 079cc8b8..16354bf4 100644 --- a/tutorials/plot_08-GroupFits.py +++ b/tutorials/plot_08-GroupFits.py @@ -1,6 +1,6 @@ """ -08: Fitting group of spectra -============================ +08: Fitting Groups of Spectra +============================= Using the group model object to run fit models across multiple power spectra. """