gwin.models.gaussian_noise module

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

class gwin.models.gaussian_noise.GaussianNoise(variable_params, data, waveform_generator, f_lower, psds=None, f_upper=None, norm=None, **kwargs)

Bases: gwin.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)\) and for noise \(\log p(d|n)\) are given by:

\[\begin{split}\log p(d|\Theta) &= -\frac{1}{2} \sum_i \left< h_i(\Theta) - d_i | h_i(\Theta) - d_i \right> \\ \log p(d|n) &= -\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 inner product is given by:

\[\left<a | b\right> = 4\Re \int_{0}^{\infty} \frac{\tilde{a}(f) \tilde{b}(f)}{S_n(f)} \mathrm{d}f,\]

where \(S_n(f)\) is the PSD in the given detector.

Note that the log prior-weighted likelihood ratio has one fewer term than the log posterior, since the \(\left<d_i|d_i\right>\) term cancels in the likelihood ratio:

\[\log \hat{\mathcal{L}} = \log p(\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. The likelihood function of the noise,

\[p(d|n) = \frac{1}{2} \sum_i \left<d_i|d_i\right>,\]

is computed on initialization and stored as the lognl attribute.

By default, the data is assumed to be equally sampled in frequency, but unequally sampled data can be supported by passing the appropriate normalization using the norm keyword argument.

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

Parameters:

• variable_params ((tuple of) string(s)) – A tuple of parameter names that will be varied.
• waveform_generator (generator class) – A generator class that creates waveforms. This must have a generate function which takes parameter values as keyword arguments, a detectors attribute which is a dictionary of detectors keyed by their names, and an epoch which specifies the start time of the generated waveform.
• 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.
• f_lower (float) – The starting frequency to use 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.
• f_upper ({None, float}) – The ending frequency to use for computing inner products. If not provided, the minimum of the largest frequency stored in the data and a given waveform will be used.
• norm ({None, float or array}) – An extra normalization weight to apply to the inner products. Can be either a float or an array. If None, 4*data.values()[0].delta_f will be used.
• **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.waveform.generator import (FDomainDetFrameGenerator,
...                                       FDomainCBCGenerator)
>>> import gwin
>>> seglen = 4
>>> sample_rate = 2048
>>> N = seglen*sample_rate/2+1
>>> fmin = 30.
>>> m1, m2, s1z, s2z, tsig, ra, dec, pol, dist = (
...     38.6, 29.3, 0., 0., 3.1, 1.37, -1.26, 2.76, 3*500.)
>>> variable_params = ['tc']
>>> generator = FDomainDetFrameGenerator(
...     FDomainCBCGenerator, 0.,
...     variable_args=variable_params, detectors=['H1', 'L1'],
...     delta_f=1./seglen, f_lower=fmin,
...     approximant='SEOBNRv2_ROM_DoubleSpin',
...     mass1=m1, mass2=m2, spin1z=s1z, spin2z=s2z,
...     ra=ra, dec=dec, polarization=pol, distance=dist)
>>> signal = generator.generate(tc=tsig)
>>> psd = pypsd.aLIGOZeroDetHighPower(N, 1./seglen, 20.)
>>> psds = {'H1': psd, 'L1': psd}
>>> model = gwin.models.GaussianNoise(
...     variable_params, signal, generator, fmin, psds=psds)

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))
278.96
>>> print('{:.2f}'.format(model.logplr))
278.96

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: 175.57+0.00j,
H1_optimal_snrsq: 351.13,
L1_cplx_loglr: 103.40+0.00j,
L1_optimal_snrsq: 206.79,
logjacobian: 0.00,
loglikelihood: 0.00,
loglr: 278.96,
logprior: 0.00

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.62

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.62

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.arange(seglen*sample_rate)/float(sample_rate)
>>> loglrs = numpy.zeros(len(times))
>>> for (ii, t) in enumerate(times):
...     model.update(tc=t)
...     loglrs[ii] = model.loglr
>>> print('tsig: {:.3f}, time of max loglr: {:.3f}'.format(
...     tsig, times[loglrs.argmax()]))
tsig: 3.100, time of max loglr: 3.100

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 = gwin.models.GaussianNoise(variable_params,
...     signal, generator, 20., psds=psds, prior=prior)
>>> model.update(tc=tsig)
>>> print('{:.2f}'.format(model.logplr))
279.88
>>> print(',\n'.join(['{}: {:.2f}'.format(s, v)
...                   for (s, v) in sorted(model.current_stats.items())]))
H1_cplx_loglr: 175.57+0.00j,
H1_optimal_snrsq: 351.13,
L1_cplx_loglr: 103.40+0.00j,
L1_optimal_snrsq: 206.79,
logjacobian: 0.00,
loglikelihood: 0.00,
loglr: 278.96,
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.

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

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

name = 'gaussian_noise'