pycbc.inference.models package

Submodules

pycbc.inference.models.analytic module

This modules provides models that have analytic solutions for the log likelihood.

class pycbc.inference.models.analytic.TestEggbox(variable_params, **kwargs)

Bases: pycbc.inference.models.base.BaseModel

The test distribution is an ‘eggbox’ function:

\[\log \mathcal{L}(\Theta) = \left[ 2+\prod_{i=1}^{n}\cos\left(\frac{\theta_{i}}{2}\right)\right]^{5}\]

The number of dimensions is set by the number of variable_params that are passed.

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied.
• **kwargs – All other keyword arguments are passed to BaseModel.

name = 'test_eggbox'

class pycbc.inference.models.analytic.TestNormal(variable_params, mean=None, cov=None, **kwargs)

Bases: pycbc.inference.models.base.BaseModel

The test distribution is an multi-variate normal distribution.

The number of dimensions is set by the number of variable_params that are passed. For details on the distribution used, see scipy.stats.multivariate_normal.

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied.
• mean (array-like, optional) – The mean values of the parameters. If None provide, will use 0 for all parameters.
• cov (array-like, optional) – The covariance matrix of the parameters. If None provided, will use unit variance for all parameters, with cross-terms set to 0.
• **kwargs – All other keyword arguments are passed to BaseModel.

Examples

Create a 2D model with zero mean and unit variance:

>>> m = TestNormal(['x', 'y'])

Set the current parameters and evaluate the log posterior:

>>> m.update(x=-0.2, y=0.1)
>>> m.logposterior
-1.8628770664093453

See the current stats that were evaluated:

>>> m.current_stats
{'logjacobian': 0.0, 'loglikelihood': -1.8628770664093453, 'logprior': 0.0}

name = 'test_normal'

class pycbc.inference.models.analytic.TestPrior(variable_params, **kwargs)

Bases: pycbc.inference.models.base.BaseModel

Uses the prior as the test distribution.

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied. Must have length 2.
• **kwargs – All other keyword arguments are passed to BaseModel.

name = 'test_prior'

class pycbc.inference.models.analytic.TestRosenbrock(variable_params, **kwargs)

Bases: pycbc.inference.models.base.BaseModel

The test distribution is the Rosenbrock function:

\[\log \mathcal{L}(\Theta) = -\sum_{i=1}^{n-1}[ (1-\theta_{i})^{2}+100(\theta_{i+1} - \theta_{i}^{2})^{2}]\]

The number of dimensions is set by the number of variable_params that are passed.

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied.
• **kwargs – All other keyword arguments are passed to BaseModel.

name = 'test_rosenbrock'

class pycbc.inference.models.analytic.TestVolcano(variable_params, **kwargs)

Bases: pycbc.inference.models.base.BaseModel

The test distribution is a two-dimensional ‘volcano’ function:

\[\Theta = \sqrt{\theta_{1}^{2} + \theta_{2}^{2}} \log \mathcal{L}(\Theta) = 25\left(e^{\frac{-\Theta}{35}} + \frac{1}{2\sqrt{2\pi}} e^{-\frac{(\Theta-5)^{2}}{8}}\right)\]

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied. Must have length 2.
• **kwargs – All other keyword arguments are passed to BaseModel.

name = 'test_volcano'

pycbc.inference.models.base module

Base class for models.

class pycbc.inference.models.base.BaseModel(variable_params, static_params=None, prior=None, sampling_transforms=None, waveform_transforms=None)

Bases: object

Base class for all models.

Given some model \(h\) with parameters \(\Theta\), Bayes Theorem states that the probability of observing parameter values \(\vartheta\) is:

\[p(\vartheta|h) = \frac{p(h|\vartheta) p(\vartheta)}{p(h)}.\]

Here:

• \(p(\vartheta|h)\) is the posterior probability;
• \(p(h|\vartheta)\) is the likelihood;
• \(p(\vartheta)\) is the prior;
• \(p(h)\) is the evidence.

This class defines properties and methods for evaluating the log likelihood, log prior, and log posteror. A set of parameter values is set using the update method. Calling the class’s log(likelihood|prior|posterior) properties will then evaluate the model at those parameter values.

Classes that inherit from this class must implement a _loglikelihood function that can be called by loglikelihood.

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied.
• static_params (dict, optional) – A dictionary of parameter names -> values to keep fixed.
• prior (callable, optional) – A callable class or function that computes the log of the prior. If None provided, will use _noprior, which returns 0 for all parameter values.
• sampling_params (list, optional) – Replace one or more of the variable_params with the given parameters for sampling.
• replace_parameters (list, optional) – The variable_params to replace with sampling parameters. Must be the same length as sampling_params.
• sampling_transforms (list, optional) – List of transforms to use to go between the variable_params and the sampling parameters. Required if sampling_params is not None.
• waveform_transforms (list, optional) – A list of transforms to convert the variable_params into something understood by the likelihood model. This is useful if the prior is more easily parameterized in parameters that are different than what the likelihood is most easily defined in. Since these are used solely for converting parameters, and not for rescaling the parameter space, a Jacobian is not required for these transforms.
• Properties –
• ---------- –
• logjacobian – Returns the log of the jacobian needed to go from the parameter space of the variable_params to the sampling params.
• logprior – Returns the log of the prior.
• loglikelihood – A function that returns the log of the likelihood function.
• logposterior – A function that returns the log of the posterior.
• loglr – A function that returns the log of the likelihood ratio.
• logplr – A function that returns the log of the prior-weighted likelihood ratio.

current_params

current_stats

Return the default_stats as a dict.

This does no computation. It only returns what has already been calculated. If a stat hasn’t been calculated, it will be returned as numpy.nan.

Returns:Dictionary of stat names -> current stat values. Return type:dict

default_stats

The stats that get_current_stats returns by default.

static extra_args_from_config(cp, section, skip_args=None, dtypes=None)

Gets any additional keyword in the given config file.

Parameters:

• cp (WorkflowConfigParser) – Config file parser to read.
• section (str) – The name of the section to read.
• skip_args (list of str, optional) – Names of arguments to skip.
• dtypes (dict, optional) – A dictionary of arguments -> data types. If an argument is found in the dict, it will be cast to the given datatype. Otherwise, the argument’s value will just be read from the config file (and thus be a string).

Returns:

Dictionary of keyword arguments read from the config file.

Return type:

dict

classmethod from_config(cp, **kwargs)

Initializes an instance of this class from the given config file.

Parameters:

• cp (WorkflowConfigParser) – Config file parser to read.
• **kwargs – All additional keyword arguments are passed to the class. Any provided keyword will over ride what is in the config file.

get_current_stats(names=None)

Return one or more of the current stats as a tuple.

This function does no computation. It only returns what has already been calculated. If a stat hasn’t been calculated, it will be returned as numpy.nan.

Parameters:names (list of str, optional) – Specify the names of the stats to retrieve. If None (the default), will return default_stats. Returns:The current values of the requested stats, as a tuple. The order of the stats is the same as the names. Return type:tuple

logjacobian

The log jacobian of the sampling transforms at the current postion.

If no sampling transforms were provided, will just return 0.

Parameters:**params – The keyword arguments should specify values for all of the variable args and all of the sampling args. Returns:The value of the jacobian. Return type:float

loglikelihood

The log likelihood at the current parameters.

This will initially try to return the current_stats.loglikelihood. If that raises an AttributeError, will call _loglikelihood` to calculate it and store it to current_stats.

logposterior

Returns the log of the posterior of the current parameter values.

The logprior is calculated first. If the logprior returns -inf (possibly indicating a non-physical point), then the loglikelihood is not called.

logprior

Returns the log prior at the current parameters.

name = None

static prior_from_config(cp, variable_params, prior_section, constraint_section)

Gets arguments and keyword arguments from a config file.

Parameters:

• cp (WorkflowConfigParser) – Config file parser to read.
• variable_params (list) – List of of model parameter names.
• prior_section (str) – Section to read prior(s) from.
• constraint_section (str) – Section to read constraint(s) from.

Returns:

The prior.

Return type:

pycbc.distributions.JointDistribution

prior_rvs(size=1, prior=None)

Returns random variates drawn from the prior.

If the sampling_params are different from the variable_params, the variates are transformed to the sampling_params parameter space before being returned.

Parameters:

• size (int, optional) – Number of random values to return for each parameter. Default is 1.
• prior (JointDistribution, optional) – Use the given prior to draw values rather than the saved prior.

Returns:

A field array of the random values.

Return type:

FieldArray

sampling_params

Returns the sampling parameters.

If sampling_transforms is None, this is the same as the variable_params.

static_params

Returns the model’s static arguments.

update(**params)

Updates the current parameter positions and resets stats.

If any sampling transforms are specified, they are applied to the params before being stored.

variable_params

Returns the model parameters.

write_metadata(fp)

Writes metadata to the given file handler.

Parameters:fp (pycbc.inference.io.BaseInferenceFile instance) – The inference file to write to.

class pycbc.inference.models.base.ModelStats

Bases: object

Class to hold model’s current stat values.

getstats(names, default=nan)

Get the requested stats as a tuple.

If a requested stat is not an attribute (implying it hasn’t been stored), then the default value is returned for that stat.

Parameters:

• names (list of str) – The names of the stats to get.
• default (float, optional) – What to return if a requested stat is not an attribute of self. Default is numpy.nan.

Returns:

A tuple of the requested stats.

Return type:

tuple

getstatsdict(names, default=nan)

Get the requested stats as a dictionary.

If a requested stat is not an attribute (implying it hasn’t been stored), then the default value is returned for that stat.

Parameters:

• names (list of str) – The names of the stats to get.
• default (float, optional) – What to return if a requested stat is not an attribute of self. Default is numpy.nan.

Returns:

A dictionary of the requested stats.

Return type:

dict

statnames

Returns the names of the stats that have been stored.

class pycbc.inference.models.base.SamplingTransforms(variable_params, sampling_params, replace_parameters, sampling_transforms)

Bases: object

Provides methods for transforming between sampling parameter space and model parameter space.

apply(samples, inverse=False)

Applies the sampling transforms to the given samples.

Parameters:

• samples (dict or FieldArray) – The samples to apply the transforms to.
• inverse (bool, optional) – Whether to apply the inverse transforms (i.e., go from the sampling args to the variable_params). Default is False.

Returns:

The transformed samples, along with the original samples.

Return type:

dict or FieldArray

classmethod from_config(cp, variable_params)

Gets sampling transforms specified in a config file.

Sampling parameters and the parameters they replace are read from the sampling_params section, if it exists. Sampling transforms are read from the sampling_transforms section(s), using transforms.read_transforms_from_config.

An AssertionError is raised if no sampling_params section exists in the config file.

Parameters:

• cp (WorkflowConfigParser) – Config file parser to read.
• variable_params (list) – List of parameter names of the original variable params.

Returns:

A sampling transforms class.

Return type:

SamplingTransforms

logjacobian(**params)

Returns the log of the jacobian needed to transform pdfs in the variable_params parameter space to the sampling_params parameter space.

Let \(\mathbf{x}\) be the set of variable parameters, \(\mathbf{y} = f(\mathbf{x})\) the set of sampling parameters, and \(p_x(\mathbf{x})\) a probability density function defined over \(\mathbf{x}\). The corresponding pdf in \(\mathbf{y}\) is then:

\[p_y(\mathbf{y}) = p_x(\mathbf{x})\left|\mathrm{det}\,\mathbf{J}_{ij}\right|,\]

where \(\mathbf{J}_{ij}\) is the Jacobian of the inverse transform \(\mathbf{x} = g(\mathbf{y})\). This has elements:

\[\mathbf{J}_{ij} = \frac{\partial g_i}{\partial{y_j}}\]

This function returns \(\log \left|\mathrm{det}\,\mathbf{J}_{ij}\right|\).

Parameters:**params – The keyword arguments should specify values for all of the variable args and all of the sampling args. Returns:The value of the jacobian. Return type:float

pycbc.inference.models.base.read_sampling_params_from_config(cp, section_group=None, section='sampling_params')

Reads sampling parameters from the given config file.

Parameters are read from the [({section_group}_){section}] section. The options should list the variable args to transform; the parameters they point to should list the parameters they are to be transformed to for sampling. If a multiple parameters are transformed together, they should be comma separated. Example:

[sampling_params]
mass1, mass2 = mchirp, logitq
spin1_a = logitspin1_a

Note that only the final sampling parameters should be listed, even if multiple intermediate transforms are needed. (In the above example, a transform is needed to go from mass1, mass2 to mchirp, q, then another one needed to go from q to logitq.) These transforms should be specified in separate sections; see transforms.read_transforms_from_config for details.

Parameters:

• cp (WorkflowConfigParser) – An open config parser to read from.
• section_group (str, optional) – Append {section_group}_ to the section name. Default is None.
• section (str, optional) – The name of the section. Default is ‘sampling_params’.

Returns:

• sampling_params (list) – The list of sampling parameters to use instead.
• replaced_params (list) – The list of variable args to replace in the sampler.

pycbc.inference.models.base_data module

Base classes for models with data.

class pycbc.inference.models.base_data.BaseDataModel(variable_params, data, recalibration=None, gates=None, **kwargs)

Bases: pycbc.inference.models.base.BaseModel

Base class for models that require data and a waveform generator.

This adds propeties for the log of the likelihood that the data contain noise, lognl, and the log likelihood ratio loglr.

Classes that inherit from this class must define _loglr and _lognl functions, in addition to the _loglikelihood requirement inherited from BaseModel.

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied.
• data (dict) – A dictionary of data, in which the keys are the detector names and the values are the data.
• recalibration (dict of pycbc.calibration.Recalibrate, optional) – Dictionary of detectors -> recalibration class instances for recalibrating data.
• gates (dict of tuples, optional) – Dictionary of detectors -> tuples of specifying gate times. The sort of thing returned by pycbc.gate.gates_from_cli.
• **kwargs – All other keyword arguments are passed to BaseModel.

data

The data that the class was initialized with.

Type:dict

lognl

Returns the log likelihood of the noise.

loglr

Returns the log of the likelihood ratio.

logplr

Returns the log of the prior-weighted likelihood ratio.

See ``BaseModel`` for additional attributes and properties.

data

Dictionary mapping detector names to data.

detectors

Returns the detectors used.

loglr

The log likelihood ratio at the current parameters.

This will initially try to return the current_stats.loglr. If that raises an AttributeError, will call _loglr` to calculate it and store it to current_stats.

lognl

The log likelihood of the model assuming the data is noise.

This will initially try to return the current_stats.lognl. If that raises an AttributeError, will call _lognl` to calculate it and store it to current_stats.

logplr

Returns the log of the prior-weighted likelihood ratio at the current parameter values.

The logprior is calculated first. If the logprior returns -inf (possibly indicating a non-physical point), then loglr is not called.

write_metadata(fp)

Adds data to the metadata that’s written.

Parameters:fp (pycbc.inference.io.BaseInferenceFile instance) – The inference file to write to.

pycbc.inference.models.gaussian_noise module

This module provides model classes that assume the noise is Gaussian.

class pycbc.inference.models.gaussian_noise.GaussianNoise(variable_params, data, low_frequency_cutoff, psds=None, high_frequency_cutoff=None, static_params=None, **kwargs)

Bases: pycbc.inference.models.base_data.BaseDataModel

Model that assumes data is stationary Gaussian noise.

With Gaussian noise the log likelihood functions for signal \(\log p(d|\Theta, h)\) and for noise \(\log p(d|n)\) are given by:

\[\begin{split}\log p(d|\Theta, h) &= \log\alpha -\frac{1}{2} \sum_i \left< d_i - h_i(\Theta) | d_i - h_i(\Theta) \right> \\ \log p(d|n) &= -\log\alpha -\frac{1}{2} \sum_i \left<d_i | d_i\right>\end{split}\]

where the sum is over the number of detectors, \(d_i\) is the data in each detector, and \(h_i(\Theta)\) is the model signal in each detector. The (discrete) inner product is given by:

\[\left<a_i | b_i\right> = 4\Re \delta f \sum_{k=k_{\mathrm{min}}}^{k_{\mathrm{max}}} \frac{\tilde{a}_i^{*}[k] \tilde{b}_i[k]}{S^(i)_n[k]},\]

where \(\delta f\) is the frequency resolution (given by 1 / the observation time \(T\)), \(k\) is an index over the discretely sampled frequencies \(f = k \delta_f\), and \(S^(i)_n[k]\) is the PSD in the given detector. The upper cutoff on the inner product \(k_{\mathrm{max}}\) is by default the Nyquist frequency k_{mathrm{max}} = N/2+1, where \(N = \lfloor T/\delta t \rfloor\) is the number of samples in the time domain, but this can be set manually to a smaller value.

The normalization factor \(\alpha\) is:

\[\alpha = \frac{1}{\left(\pi T\right)^{N/2} \prod_{k=k_\mathrm{min}}^{k_{\mathrm{max}}} S_n[k]}.\]

Note that the log likelihood ratio has fewer terms than the log likelihood, since the normalization and \(\left<d_i|d_i\right>\) terms cancel:

\[\log \mathcal{L}(\Theta) = \sum_i \left[ \left<h_i(\Theta)|d_i\right> - \frac{1}{2} \left<h_i(\Theta)|h_i(\Theta)\right> \right]\]

Upon initialization, the data is whitened using the given PSDs. If no PSDs are given the data and waveforms returned by the waveform generator are assumed to be whitened.

For more details on initialization parameters and definition of terms, see models.BaseDataModel.

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied.
• data (dict) – A dictionary of data, in which the keys are the detector names and the values are the data (assumed to be unwhitened). The list of keys must match the waveform generator’s detectors keys, and the epoch of every data set must be the same as the waveform generator’s epoch.
• low_frequency_cutoff (dict) – A dictionary of starting frequencies, in which the keys are the detector names and the values are the starting frequencies for the respective detectors to be used for computing inner products.
• psds ({None, dict}) – A dictionary of FrequencySeries keyed by the detector names. The dictionary must have a psd for each detector specified in the data dictionary. If provided, the inner products in each detector will be weighted by 1/psd of that detector.
• high_frequency_cutoff (dict, optional) – A dictionary of ending frequencies, in which the keys are the detector names and the values are the ending frequencies for the respective detectors to be used for computing inner products. If not provided, the minimum of the largest frequency stored in the data and a given waveform will be used.
• static_params (dict, optional) – A dictionary of parameter names -> values to keep fixed.
• **kwargs – All other keyword arguments are passed to BaseDataModel.

Examples

Create a signal, and set up the model using that signal:

>>> from pycbc import psd as pypsd
>>> from pycbc.inference.models import GaussianNoise
>>> from pycbc.waveform.generator import (FDomainDetFrameGenerator,
...                                       FDomainCBCGenerator)
>>> seglen = 4
>>> sample_rate = 2048
>>> N = seglen*sample_rate/2+1
>>> fmin = 30.
>>> static_params = {'approximant': 'IMRPhenomD', 'f_lower': fmin,
...                  'mass1': 38.6, 'mass2': 29.3,
...                  'spin1z': 0., 'spin2z': 0., 'ra': 1.37, 'dec': -1.26,
...                  'polarization': 2.76, 'distance': 3*500.}
>>> variable_params = ['tc']
>>> tsig = 3.1
>>> generator = FDomainDetFrameGenerator(
...     FDomainCBCGenerator, 0., detectors=['H1', 'L1'],
...     variable_args=variable_params,
...     delta_f=1./seglen, **static_params)
>>> signal = generator.generate(tc=tsig)
>>> psd = pypsd.aLIGOZeroDetHighPower(N, 1./seglen, 20.)
>>> psds = {'H1': psd, 'L1': psd}
>>> low_frequency_cutoff = {'H1': fmin, 'L1': fmin}
>>> model = GaussianNoise(variable_params, signal, low_frequency_cutoff,
                          psds=psds, static_params=static_params)

Set the current position to the coalescence time of the signal:

>>> model.update(tc=tsig)

Now compute the log likelihood ratio and prior-weighted likelihood ratio; since we have not provided a prior, these should be equal to each other:

>>> print('{:.2f}'.format(model.loglr))
282.43
>>> print('{:.2f}'.format(model.logplr))
282.43

Print all of the default_stats:

>>> print(',\n'.join(['{}: {:.2f}'.format(s, v)
...                   for (s, v) in sorted(model.current_stats.items())]))
H1_cplx_loglr: 177.76+0.00j,
H1_optimal_snrsq: 355.52,
L1_cplx_loglr: 104.67+0.00j,
L1_optimal_snrsq: 209.35,
logjacobian: nan,
loglikelihood: 835680.31,
loglr: 282.43,
logprior: nan

Compute the SNR; for this system and PSD, this should be approximately 24:

>>> from pycbc.conversions import snr_from_loglr
>>> x = snr_from_loglr(model.loglr)
>>> print('{:.2f}'.format(x))
23.77

Since there is no noise, the SNR should be the same as the quadrature sum of the optimal SNRs in each detector:

>>> x = (model.det_optimal_snrsq('H1') +
...      model.det_optimal_snrsq('L1'))**0.5
>>> print('{:.2f}'.format(x))
23.77

Using the same model, evaluate the log likelihood ratio at several points in time and check that the max is at tsig:

>>> import numpy
>>> times = numpy.linspace(tsig-1, tsig+1, num=101)
>>> loglrs = numpy.zeros(len(times))
>>> for (ii, t) in enumerate(times):
...     model.update(tc=t)
...     loglrs[ii] = model.loglr
>>> print('tsig: {:.2f}, time of max loglr: {:.2f}'.format(
...     tsig, times[loglrs.argmax()]))
tsig: 3.10, time of max loglr: 3.10

Create a prior and use it (see distributions module for more details):

>>> from pycbc import distributions
>>> uniform_prior = distributions.Uniform(tc=(tsig-0.2,tsig+0.2))
>>> prior = distributions.JointDistribution(variable_params, uniform_prior)
>>> model = pycbc.inference.models.GaussianNoise(variable_params,
...     signal, low_frequency_cutoff, psds=psds, prior=prior,
...     static_params=static_params)
>>> model.update(tc=tsig)
>>> print('{:.2f}'.format(model.logplr))
283.35
>>> print(',\n'.join(['{}: {:.2f}'.format(s, v)
...                   for (s, v) in sorted(model.current_stats.items())]))
H1_cplx_loglr: 177.76+0.00j,
H1_optimal_snrsq: 355.52,
L1_cplx_loglr: 104.67+0.00j,
L1_optimal_snrsq: 209.35,
logjacobian: 0.00,
loglikelihood: 835680.31,
loglr: 282.43,
logprior: 0.92

det_cplx_loglr(det)

Returns the complex log likelihood ratio in the given detector.

Parameters:det (str) – The name of the detector. Returns:The complex log likelihood ratio. Return type:complex float

det_lognl(det)

Returns the log likelihood of the noise in the given detector:

\[\log p(d_i|n_i) = \log \alpha_i - \frac{1}{2} \left<d_i | d_i\right>.\]

Parameters:det (str) – The name of the detector. Returns:The log likelihood of the noise in the requested detector. Return type:float

det_lognorm(det)

The log of the likelihood normalization in the given detector.

det_optimal_snrsq(det)

Returns the opitmal SNR squared in the given detector.

Parameters:det (str) – The name of the detector. Returns:The opimtal SNR squared. Return type:float

classmethod from_config(cp, **kwargs)

Initializes an instance of this class from the given config file.

Parameters:

• cp (WorkflowConfigParser) – Config file parser to read.
• **kwargs – All additional keyword arguments are passed to the class. Any provided keyword will over ride what is in the config file.

high_frequency_cutoff

The high frequency cutoff of the inner product.

lognorm

The log of the normalization of the log likelihood.

low_frequency_cutoff

The low frequency cutoff of the inner product.

name = 'gaussian_noise'

psds

Returns the psds that are set.

waveform_generator

The waveform generator used.

write_metadata(fp)

Adds writing the psds and lognl, since it’s a constant.

The lognl is written to the sample group’s attrs.

Parameters:fp (pycbc.inference.io.BaseInferenceFile instance) – The inference file to write to.

pycbc.inference.models.gaussian_noise.create_waveform_generator(variable_params, data, recalibration=None, gates=None, **static_params)

Creates a waveform generator for use with a model.

Parameters:

• variable_params (list of str) – The names of the parameters varied.
• data (dict) – Dictionary mapping detector names to either a <pycbc.types.TimeSeries TimeSeries> or <pycbc.types.FrequencySeries FrequencySeries>.
• recalibration (dict, optional) – Dictionary mapping detector names to <pycbc.calibration.Recalibrate> instances for recalibrating data.
• gates (dict of tuples, optional) – Dictionary of detectors -> tuples of specifying gate times. The sort of thing returned by pycbc.gate.gates_from_cli().

Returns:

A waveform generator for frequency domain generation.

Return type:

pycbc.waveform.FDomainDetFrameGenerator

pycbc.inference.models.gaussian_noise.high_frequency_cutoff_from_config(cp)

Gets the high frequency cutoff, if one is provided for a detector from the given config file for likelihood computation.

This looks for options IFO-high-frequency-cutoff in the [model] section, where IFO is the name of the detector, and casts it to float.

Parameters:cp (WorkflowConfigParser) – Config file parser to read. Returns:Dictionary with the IFO names as the keys and the respective high frequency cutoffs to be used in the likelihood calculation as the values. If high frequency cutoffs for no detectors are provided, returns an empty dictionary. Return type:dict or None

pycbc.inference.models.gaussian_noise.low_frequency_cutoff_from_config(cp)

Gets the low frequency cutoff for all detectors to be used from the given config file.

The low-frequency-cutoff for each detector should be provided using an option IFO-low-frequency-cutoff in the [model] section, where IFO is the name of the detector. The low frequency cutoff value is then casted to float. If the casting to float fails, an error is raised.

Parameters:cp (WorkflowConfigParser) – Config file parser to read. Returns:Dictionary with the IFO names as the keys and the respective low frequency cutoffs to be used in the likelihood calculation as the values. Return type:dict

pycbc.inference.models.marginalized_gaussian_noise module

This module provides model classes that assume the noise is Gaussian and allows for the likelihood to be marginalized over phase and/or time and/or distance.

class pycbc.inference.models.marginalized_gaussian_noise.MarginalizedGaussianNoise(variable_params, data, waveform_generator, f_lower, psds=None, f_upper=None, norm=None, time_marginalization=False, distance_marginalization=False, phase_marginalization=False, marg_prior=None, **kwargs)

Bases: pycbc.inference.models.gaussian_noise.GaussianNoise

The likelihood is analytically marginalized over phase and/or time and/or distance.

For the case of marginalizing over phase, the signal, can be written as:

\[\tilde{h}(f; \Theta, \phi) = A(f; \Theta)e^{i\Psi(f; \Theta) + i \phi},\]

where \(\phi\) is an arbitrary phase constant. This phase constant can be analytically marginalized over with a uniform prior as follows: assuming the noise is stationary and Gaussian (see GaussianNoise for details), the posterior is:

\[\begin{split}p(\Theta,\phi|d) &\propto p(\Theta)p(\phi)p(d|\Theta,\phi) \\ &\propto p(\Theta)\frac{1}{2\pi}\exp\left[ -\frac{1}{2}\sum_{i}^{N_D} \left< h_i(\Theta,\phi) - d_i, h_i(\Theta,\phi) - d_i \right>\right].\end{split}\]

Here, the sum is over the number of detectors \(N_D\), \(d_i\) and \(h_i\) are the data and signal in detector \(i\), respectively, and we have assumed a uniform prior on \(phi \in [0, 2\pi)\). With the form of the signal model given above, the inner product in the exponent can be written as:

\[\begin{split}-\frac{1}{2}\left<h_i - d_i, h_i- d_i\right> &= \left<h_i, d_i\right> - \frac{1}{2}\left<h_i, h_i\right> - \frac{1}{2}\left<d_i, d_i\right> \\ &= \Re\left\{O(h^0_i, d_i)e^{-i\phi}\right\} - \frac{1}{2}\left<h^0_i, h^0_i\right> - \frac{1}{2}\left<d_i, d_i\right>,\end{split}\]

where:

\[\begin{split}h_i^0 &\equiv \tilde{h}_i(f; \Theta, \phi=0); \\ O(h^0_i, d_i) &\equiv 4 \int_0^\infty \frac{\tilde{h}_i^*(f; \Theta,0)\tilde{d}_i(f)}{S_n(f)}\mathrm{d}f.\end{split}\]

Gathering all of the terms that are not dependent on \(\phi\) together:

\[\alpha(\Theta, d) \equiv \exp\left[-\frac{1}{2}\sum_i \left<h^0_i, h^0_i\right> + \left<d_i, d_i\right>\right],\]

we can marginalize the posterior over \(\phi\):

\[\begin{split}p(\Theta|d) &\propto p(\Theta)\alpha(\Theta,d)\frac{1}{2\pi} \int_{0}^{2\pi}\exp\left[\Re \left\{ e^{-i\phi} \sum_i O(h^0_i, d_i) \right\}\right]\mathrm{d}\phi \\ &\propto p(\Theta)\alpha(\Theta, d)\frac{1}{2\pi} \int_{0}^{2\pi}\exp\left[ x(\Theta,d)\cos(\phi) + y(\Theta, d)\sin(\phi) \right]\mathrm{d}\phi.\end{split}\]

The integral in the last line is equal to \(2\pi I_0(\sqrt{x^2+y^2})\), where \(I_0\) is the modified Bessel function of the first kind. Thus the marginalized log posterior is:

\[\log p(\Theta|d) \propto \log p(\Theta) + I_0\left(\left|\sum_i O(h^0_i, d_i)\right|\right) - \frac{1}{2}\sum_i\left[ \left<h^0_i, h^0_i\right> - \left<d_i, d_i\right> \right]\]

For the case of marginalizing over distance, the signal can be written as,

\[\tilde{h}_{j} = \frac{1}{D} \tilde{h}_{j}^{0}\]

The distance can be analytically marginalized over with a uniform prior as follows: assuming the noise is stationary and Gaussian (see GaussianNoise for details), the likelihood is:

\[log L = -\frac{1}{2}\left<d-h|d-h\right>\]

We see that :math: <h|h> is inversely proportional to distance squared and :math: <h|d> is inversely proportional to distance. The log likelihood is therefore

\[\log L = -\frac{1}{2}\left<d|d\right> - \frac{1}{2D^{2}} \left<h|h\right> + \frac{1}{D}\left<h|d\right>\]

Consequently, the likelihood marginalised over distance is simply

\[\log L = \log\left(\int_{0}^{D}{L p(D) dD}\right)\]

If we assume a flat prior

\[\log L = \log\left(\int_{0}^{D}{\exp{\log L} dD}\right)\]

For the case of marginalizing over time, the signal can be written as,

\[\tilde{h}_{j} = \tilde{h}_{j}^{0} \exp\left(-2\pi ij\Delta ft\right)\]

The time can be analytically marginalized over with a uniform prior as follows: assuming the noise is stationary and Gaussian (see GaussianNoise for details), the likelihood is:

\[\log L = -\frac{1}{2}\left<d-h|d-h\right>\]

We note that :math: <h|h> and :math: <d|d> are time independent while :math: <d|h> is dependent of time

\[\left<d|h\right>(t) = 4\Delta f\sum_{j=0}^{N/2} \frac{\tilde{d}_{j}^{*} \tilde{h}_{j}^{0}}{S_{j}} \exp\left(-2\pi ij \Delta f t\right)\]

For integer timesteps :math: t=kDelta t

\[ \begin{align}\begin{aligned}\left<d|h\right>(k\Delta t) = 4\Delta f\sum_{j=0}^{N/2} \frac{ \tilde{d}_{j}^{*}\tilde{h}_{j}^{0}} {S_{j}}\exp(-2\pi \frac{ijk}{N}\\\left<d|h\right>(k\Delta t) = 2\Delta f\sum_{j=0}^{N} \frac{ \tilde{d}_{j}^{*}\tilde{h}_{j}^{0}} {S_{j}} \exp(-2\pi \frac{ijk}{N}\end{aligned}\end{align} \]

Using a FFT, this expression can be evaluated efficiently for all :math: k

\[\left<d|h\right>(k\Delta t) = 2\Delta f FFT_{k} (\frac{dh}{S})\]

since :math:: left<h|dright> = left<d|hright>^{*},

\[\left<d|h\right> + \left<h|d\right> = 4\Delta f FFT_{k} (\frac{dh}{S})\]

and so the likelihood marginalised over time is simply

\[\log{L} = \log\left(\int_{0}^{T} np.exp(\np.log(L) p(t))\right)\]

where p(t) is the prior. If we assume a flat prior then,

\[\log{L} = \log\left(\int_{0}^{T} np.exp(\np.log(L))\right)\]

Parameters:

• time_marginalization (bool, optional) – A Boolean operator which determines if the likelihood is marginalized over time
• phase_marginalization (bool, optional) – A Boolean operator which determines if the likelihood is marginalized over phase
• distance_marginalization (bool, optional) – A Boolean operator which determines if the likelihood is marginalized over distance
• marg_prior (list, optional) – An instance of pycbc.distributions which returns a list of prior distributions to be used when marginalizing the likelihood
• **kwargs – All other keyword arguments are passed to GaussianNoise.

classmethod from_config(cp, data=None, delta_f=None, delta_t=None, gates=None, recalibration=None, **kwargs)

Initializes an instance of this class from the given config file.

Parameters:

• cp (WorkflowConfigParser) – Config file parser to read.
• data (dict) – A dictionary of data, in which the keys are the detector names and the values are the data. This is not retrieved from the config file, and so must be provided.
• delta_f (float) – The frequency spacing of the data; needed for waveform generation.
• delta_t (float) – The time spacing of the data; needed for time-domain waveform generators.
• recalibration (dict of pycbc.calibration.Recalibrate, optional) – Dictionary of detectors -> recalibration class instances for recalibrating data.
• gates (dict of tuples, optional) – Dictionary of detectors -> tuples of specifying gate times. The sort of thing returned by pycbc.gate.gates_from_cli.
• **kwargs – All additional keyword arguments are passed to the class. Any provided keyword will over ride what is in the config file.

name = 'marginalized_gaussian_noise'

class pycbc.inference.models.marginalized_gaussian_noise.MarginalizedPhaseGaussianNoise(variable_params, data, low_frequency_cutoff, psds=None, high_frequency_cutoff=None, static_params=None, **kwargs)

Bases: pycbc.inference.models.gaussian_noise.GaussianNoise

The likelihood is analytically marginalized over phase.

This class can be used with signal models that can be written as:

\[\tilde{h}(f; \Theta, \phi) = A(f; \Theta)e^{i\Psi(f; \Theta) + i \phi},\]

where \(\phi\) is an arbitrary phase constant. This phase constant can be analytically marginalized over with a uniform prior as follows: assuming the noise is stationary and Gaussian (see GaussianNoise for details), the posterior is:

\[\begin{split}p(\Theta,\phi|d) &\propto p(\Theta)p(\phi)p(d|\Theta,\phi) \\ &\propto p(\Theta)\frac{1}{2\pi}\exp\left[ -\frac{1}{2}\sum_{i}^{N_D} \left< h_i(\Theta,\phi) - d_i, h_i(\Theta,\phi) - d_i \right>\right].\end{split}\]

Here, the sum is over the number of detectors \(N_D\), \(d_i\) and \(h_i\) are the data and signal in detector \(i\), respectively, and we have assumed a uniform prior on \(phi \in [0, 2\pi)\). With the form of the signal model given above, the inner product in the exponent can be written as:

\[\begin{split}-\frac{1}{2}\left<h_i - d_i, h_i- d_i\right> &= \left<h_i, d_i\right> - \frac{1}{2}\left<h_i, h_i\right> - \frac{1}{2}\left<d_i, d_i\right> \\ &= \Re\left\{O(h^0_i, d_i)e^{-i\phi}\right\} - \frac{1}{2}\left<h^0_i, h^0_i\right> - \frac{1}{2}\left<d_i, d_i\right>,\end{split}\]

where:

\[\begin{split}h_i^0 &\equiv \tilde{h}_i(f; \Theta, \phi=0); \\ O(h^0_i, d_i) &\equiv 4 \int_0^\infty \frac{\tilde{h}_i^*(f; \Theta,0)\tilde{d}_i(f)}{S_n(f)}\mathrm{d}f.\end{split}\]

Gathering all of the terms that are not dependent on \(\phi\) together:

\[\alpha(\Theta, d) \equiv \exp\left[-\frac{1}{2}\sum_i \left<h^0_i, h^0_i\right> + \left<d_i, d_i\right>\right],\]

we can marginalize the posterior over \(\phi\):

\[\begin{split}p(\Theta|d) &\propto p(\Theta)\alpha(\Theta,d)\frac{1}{2\pi} \int_{0}^{2\pi}\exp\left[\Re \left\{ e^{-i\phi} \sum_i O(h^0_i, d_i) \right\}\right]\mathrm{d}\phi \\ &\propto p(\Theta)\alpha(\Theta, d)\frac{1}{2\pi} \int_{0}^{2\pi}\exp\left[ x(\Theta,d)\cos(\phi) + y(\Theta, d)\sin(\phi) \right]\mathrm{d}\phi.\end{split}\]

The integral in the last line is equal to \(2\pi I_0(\sqrt{x^2+y^2})\), where \(I_0\) is the modified Bessel function of the first kind. Thus the marginalized log posterior is:

\[\log p(\Theta|d) \propto \log p(\Theta) + I_0\left(\left|\sum_i O(h^0_i, d_i)\right|\right) - \frac{1}{2}\sum_i\left[ \left<h^0_i, h^0_i\right> - \left<d_i, d_i\right> \right]\]

name = 'marginalized_phase'

pycbc.inference.models.single_template module

This module provides model classes that assume the noise is Gaussian.

class pycbc.inference.models.single_template.SingleTemplate(data, psds, low_frequency_cutoff=None, high_frequency_cutoff=None, sample_rate=32768, **kwargs)

Bases: pycbc.inference.models.base_data.BaseDataModel

Model that assumes we know all the intrinsic parameters.

This model assumes we know all the intrinsic parameters, and are only maximizing over the extrinsic ones. We also assume a dominant mode waveform approximant only and non-precessing.

name = 'single_template'

Module contents

This package provides classes and functions for evaluating Bayesian statistics assuming various noise models.

class pycbc.inference.models.CallModel(model, callstat, return_all_stats=True)

Bases: object

Wrapper class for calling models from a sampler.

This class can be called like a function, with the parameter values to evaluate provided as a list in the same order as the model’s variable_params. In that case, the model is updated with the provided parameters and then the callstat retrieved. If return_all_stats is set to True, then all of the stats specified by the model’s default_stats will be returned as a tuple, in addition to the stat value.

The model’s attributes are promoted to this class’s namespace, so that any attribute and method of model may be called directly from this class.

This class must be initalized prior to the creation of a Pool object.

Parameters:

• model (Model instance) – The model to call.
• callstat (str) – The statistic to call.
• return_all_stats (bool, optional) – Whether or not to return all of the other statistics along with the callstat value.

Examples

Create a wrapper around an instance of the TestNormal model, with the callstat set to logposterior:

>>> from pycbc.inference.models import TestNormal, CallModel
>>> model = TestNormal(['x', 'y'])
>>> call_model = CallModel(model, 'logposterior')

Now call on a set of parameter values:

>>> call_model([0.1, -0.2])
(-1.8628770664093453, (0.0, 0.0, -1.8628770664093453))

Note that a tuple of all of the model’s default_stats were returned in addition to the logposterior value. We can shut this off by toggling return_all_stats:

>>> call_model.return_all_stats = False
>>> call_model([0.1, -0.2])
-1.8628770664093453

Attributes of the model can be called from the call model. For example:

>>> call_model.variable_params
('x', 'y')

pycbc.inference.models.read_from_config(cp, **kwargs)

Initializes a model from the given config file.

The section must have a name argument. The name argument corresponds to the name of the class to initialize.

Parameters:

• cp (WorkflowConfigParser) – Config file parser to read.
• **kwargs – All other keyword arguments are passed to the from_config method of the class specified by the name argument.

Returns:

The initialized model.

Return type:

cls