PgvLikelihood

class kszx.PgvLikelihood(data, surrs, bias_matrix, jeffreys_prior=False)

Likelihood function for one or more kSZ power spectra of the form \(P_{gv}(k)\).

This class can be used to fit \((f_{NL}, b_v)\) values or run MCMCs. Currently limited to using \(P_{gv}(k)\) only – it would be interesting to extend the analysis to include \(P_{gg}(k)\) and/or \(P_{vv}(k)\).

The PgvLikelihood constructor has an abstract syntax, which may not be the most convenient. Instead of calling the constructor directly, you may want to call from_pipeline_outdir() as follows:

dirname = ...   # directory name containing pipeline outputs
pout = kszx.KszPipeOutdir(dirname)

# The 'fields=[[1,0]]' argument selects 90 GHz.
# See the from_pipeline_outdir() docstring for more examples.
lik = kszx.PgvLikelihood.from_pipeline_outdir(pout, fields=[[1,0]], nkbins=10)

The rest of this docstring describes the PgvLikelihood constructor syntax (even though, as explained above, you probably don’t want to call the constructor directly!) First we define:

  • B = number of bias parameters in the MCMC (usually 1)

  • V = number of velocity reconstructions (usually 1)

  • K = number of k-bins

  • S = number of surrogate sims

  • D = V*K = number of components in “data vector”.

Note the distinction between a 90+150 GHz “sum” fit (V=1), and a 90+150 GHz “joint” fit (V=2). In the first case, we add the 90+150 bandpowers and construct a likelihood based on their sum, and in the second case we do a joint fit to both sets of bandpowers (i.e. doubling the size of the data vector and covariance matrices).

Then the constructor arguments are as follows:

  • data: \(P_{gv}(k)\) “data” bandpowers, an array of shape (V,K).

  • surrs: \(P_{gv}(k)\) surrogate bandpowers, an array of shape (S,2,2,V,K). The first length-2 index is an fnl exponent in {0,1}. The second length-2 index is a bv exponent in {0,1}.

  • bias_matrix: array of shape \((B,V)\), where \(B \le V\). This gives the correspondence between bias parameters and velocity reconstruction. There are basically 3 cases of interest here:

    • bias_matrix = [[1,1]]: analysis of single velocity reconstruction field (V=1).

    • bias_matrix = [[1,1]]: joint analysis of two velocity reconstruction fields (V=2), with bias assumed to be the same for both freqs.

    • bias_matrix = [[1,0],[0,1]]: joint analysis of two velocity reconstruction fields (V=2), with two independent bias parameters.

static from_pipeline_outdir(pout, fields, nkbins, multi_bias=None, jeffreys_prior=None)

Constructs a PgvLikelhood from a KszPipeOutdir instance.

Usually more convenient than calling the PgvLikelihood constructor directly. Example usage:

dirname = ...   # directory name containing pipeline outputs
pout = kszx.KszPipeOutdir(dirname)

# The 'fields=[[1,0]]' argument selects 90 GHz.
# See below for more examples.
lik = kszx.PgvLikelihood.from_pipeline_outdir(pout, fields=[[1,0]], nkbins=10)

The pout argument is a KszPipeOutdir object.

The fields argument is a V-by-2 matrix. Each row of the matrix selects a linear combination of the 90+150 GHz velocity reconstructions. For example:

  • fields = [[1,0]] for 90 GHz analysis

  • fields = [[0,1]] for 150 GHz analysis

  • fields = [[1,1]] for (90+150) “sum map” analysis

  • fields = [[1,-1]] for null (90-150) “null map” analysis

  • fields = [[1,0],[0,1]] for joint analysis with both (90+150 GHz) sets of bandpowers

The multi_bias argument is only needed for a joint analysis (i.e. \(V > 1\)) and determines whether each freq channel has an independent bias parameter (multi_bias=True), or whether all frequency channels have the same bias (multi_bias=False).

specialize_surrogates(fnl, bv, flatten)

Returns a shape (S,V,K) array, by “specializing” surrogates to specified \((f_{NL}, b_v)\).

The bv arugment should be an array of shape (B,) where \(B\) is the number of bias parameters in the likelihood. If \(B=1\), then bv can be a scalar.

Convenient but slow! Used in many places, but not fast_likelihood().

slow_mean_and_cov(fnl, bv)

Computes mean and covaraince of surrogates, for specified \((f_{NL}, b_v)\).

Returns (mean, cov) where mean.shape=(D,) and cov.shape=(D,D).

The bv arugment should be an array of shape (B,) where \(B\) is the number of bias parameters in the likelihood. If \(B=1\), then bv can be a scalar.

Brute force implementation: we call kszx.PgvLikelihood.specialize_surrogates()`(), then compute the mean and covariance with np.mean() and np.cov().

slow_mean_and_cov_gradients(fnl, bv)

Computes the gradient of the mean and covariance with respect to \((f_{NL}, b_v)\).

Returns grad_mean, grad_cov, where both gradients are represented as arrays with an extra length-(B+1) axis, i.e.

  • mean.shape = (D,)  =>  grad_mean.shape = (B+1,D)

  • cov.shape = (D,D)  =>  grad_cov.shape = (B+1,D,D)

Uses boneheaded algorithm: since mean, cov are at most quadratic in \(f_{NL}\) and \(b_v\), naive finite difference is exact (and independent of step sizes).

The bv arugment should be an array of shape (B,) where \(B\) is the number of bias parameters in the likelihood. If \(B=1\), then bv can be a scalar.

This method is only used to compute log-likelihoods with the Jeffreys prior.

slow_log_likelihood(fnl, bv)

Returns the log-likelihood at the specified \((f_{NL}, b_v)\).

The bv arugment should be an array of shape (B,) where \(B\) is the number of bias parameters in the likelihood. If \(B=1\), then bv can be a scalar.

run_mcmc(nwalkers=8, nsamples=10000, discard=1000, thin=5)

Initializes self.samples to an array of shape (N,B+1) where N is large.

show_mcmc(title=None)

Makes a corner plot from MCMC results. Intended to be called from jupyter.

analyze_chi2(fnl, bv, ddof=None)

Computes a \(\chi^2\) statistic, which compares \(P_{gv}^{data}(k)\) to model with given \((f_{NL}, b_v)\).

Returns (\(\chi^2\), \(N_{dof}\), \(p\)-value).

The bv arugment should be an array of shape (B,) where \(B\) is the number of bias parameters in the likelihood. If \(B=1\), then bv can be a scalar.

The ddof argument is used to compute the number of degrees of freedom, as ndof = nkbins - ddof. If ddof=None, then it will be equal to the number of nonzero (fnl, bias) params. (This is usually the correct choice.)

fit_fnl_and_bv(fnl0=0, bv0=0.3)

Returns \((f_{NL}, b_v)\) obtained by maximizing joint likelihood.

Note that \(b_v\) is represented as a 1-d array with shape (B,), even if \(B=1\).

fit_bv(fnl=0, bv0=0.3)

Returns \(b_v\) obtained by maximizing conditional likelihood at the given $f_{NL}.

Note that \(b_v\) is represented as a 1-d array with shape (B,), even if \(B=1\).

compute_snr()

Returns total SNR of \(P_{gv}(k)\). Does not assume a model for \(P_{gv}(k)\).

fast_mean_and_cov(fnl, bv, grad=False)

Equivalent to slow_mean_and_cov(), but faster. Intended for use in MCMC.

Note that bv must be a 1-d array of length B, i.e. scalar is not allowed.

fast_log_likelihood(fnl, bv)

Equivalent to slow_log_likelihood(), but faster. Intended for use in MCMC.

Note that bv must be a 1-d array of length B, i.e. scalar is not allowed.

test_fast_mean_and_cov()

Test fast_mean_and_cov(), by checking that it agrees with slow_mean_and_cov() at 10 random points.

test_fast_likelihood()

Test fast_log_likelihood(), by checking that it agrees with slow_log_likelihood() at 10 random points.

static make_random()

Construct and return a PgvLikelihood with random (data, surrs, bias_matrix).

Useful for standalone testing of ‘class PgvLikelihood’, in order to construct an “interesting” PgvLikelihood, in a situation where KSZ pipeline outputs are not available.

static run_tests()

Runs standalone tests of ‘class PgvLikelihood’. (Where “standalone” means that no KSZ pipeline outputs are needed.)