peakbagging

The peakbagging module contains the classes used for the peakbag stage of PBjam. The main class to interact with is the peakbag class, which will handle potential slicing of the spectrum or not and combine results for each case.

class pbjam.peakbagging.DynestyPeakbag(*args, **kwargs)

class pbjam.peakbagging.EmceePeakbag(*args, **kwargs)

class pbjam.peakbagging.basePeakbag(f, s, ell, freq, height, width, zeta, dnu, d02, freqLimits, rotAsym, RV, addPriors={}, sampling='emcee', **kwargs)

Base model class for peakbagging a section of the power spectrum.

Parameters:

• f (array-like) – The frequency array of the spectrum.

• s (array-like) – The values of the power density spectrum.

• ell (array-like) – Angular degree of the modes.

• freq (array-like) – Frequencies of the modes corresponding to the angular degrees.

• height (array-like) – Heights of the modes corresponding to the angular degrees. Default is None, in which case an SNR of 1 is assumed. Providing betters estimates from modeID is however strongly recommended.

• width (array-like) – Widths of the modes corresponding to the angular degrees. Default is None, in which case an width of 0.1 is assumed. Providing betters estimates from modeID is however strongly recommended.

• zeta (array-like) – Mixed-mode coupling factors of the modes corresponding to the angular degrees.

• dnu (float) – Large frequency separation.

• d02 (float) – Small frequency separation between l=0 and l=2 modes.

• freqLimits (list) – Frequency limits for the peakbagging.

• rotAsym (array-like) – Rotational asymmetry parameters.

• RV (array-like) – Radial velocity of the star in km/s.

• addPriors (dict, optional) – Additional priors to be added. Default is an empty dictionary.

• **kwargs (dict) – Additional keyword arguments.

AddLikeTerms(theta, thetaU)

Add the additional probabilities to likelihood

Adds the additional observational data likelihoods to the PSD likelihood.

Parameters:

p (list) – Sampling parameters.

Returns:

lnp – The likelihood of a sample given the parameter PDFs.

Return type:

float

chi_sqr(mod)

Chi^2 2 dof likelihood

Evaluates the likelihood of observing the data given the model.

Parameters:

mod (jax device array) – Spectrum model.

Returns:

L – Likelihood of the data given the model

Return type:

float

lnlikelihood(theta)

Likelihood function for set of model parameters

Evaluates the likelihood function for a set of model parameters given the data. This includes the constraint from the observed variables.

The samples l are drawn from the latent parameter priors and are first projected into the model space before the model is constructed and the likelihood is constructed.

Parameters:

l (list) – Array of latent parameters

Returns:

lnlike – The log likelihood evaluated at the model parameters p.

Return type:

float

model(thetaU)

Computes the model spectrum using a sum of Lorentzian profiles with a multiplet of 2l+1 modes per mode nl.

The rotational splitting is a weighted sum of the core and envelope rotation rates, where the weighting is the degree of mixing zeta. The splitting is symmetric.

Parameters:

thetaU (dict) – A dictionary of model parameters.

Returns:

modes – The computed model spectrum.

Return type:

array-like

parseSamples(samples, N=10000)

Parses the samples to extract and organize the model parameters.

Attempts to include at most N samples from the model, but will default to the actual number of samples of the model parameters if it’s less than N.

The resulting dictionary contains some global parameters, ell, enn, emm etc. and two dictionaries, one containing the samples drawn and one with their summary statistics.

Parameters:

• smp (dict) – A dictionary of sampled parameters.

• Nmax (int, optional) – Maximum number of samples to include. Default is 5000.

Returns:

result – A dictionary containing parsed and organized model parameters.

Return type:

dict

setAddLikeTerms()

Set attribute containing additional observational data

Additional observational data other than the power spectrum goes here.

Can be Teff or bp_rp color, but may also be additional constraints on e.g., numax, dnu.

setFreqRange()

Get frequency range around numax for model

Returns:

idx – Array of boolean values defining the interval of the frequency axis where the oscillation modes present.

Return type:

jax device array

setLabels()

Sets labels for the model parameters.

For each key in self.variables, if the key is ‘freq’, ‘height’, or ‘width’, it appends numbered labels for each mode.

setPriors(freq_err=0.03)

Sets the prior distributions for model parameters.

Parameters:

freq_err (float, optional) – The error in frequency, used to set the scale in percent of dnu of the normal distribution for frequency priors. Default is 0.03.

Notes

• Initializes the priors dictionary and updates it with additional priors specified in self.addPriors.

• For each mode, sets normal priors for ‘freq’, ‘height’, and ‘width’ if not already specified in self.priors.

• Sets uniform priors for ‘nurot_e’ (envelope rotation) and ‘nurot_c’ (core rotation) if not already specified.

• Sets a truncated sine prior for ‘inc’ (inclination) if not already specified.

• Sets a normal prior for ‘shot’ noise if not already specified.

Raises:

ValueError if the length of labels does not match the length of priors. –

unpackParams(theta)

Cast the parameters in a dictionary

Parameters:

theta (array) – Array of parameters drawn from the posterior distribution.

Returns:

thetaU – The unpacked parameters.

Return type:

dict

unpackSamples(samples)

Unpack a set of parameter samples into a dictionary of arrays.

Parameters:

samples (array-like) – A 2D array of shape (n, m), where n is the number of samples and m is the number of parameters.

Returns:

S – A dictionary containing the parameter values for each parameter label.

Return type:

dict

class pbjam.peakbagging.jointRotInc(pb, NKDE=2500, bw=0.03)

A class to perform joint posterior sampling for rotation and inclination parameters using emcee.

Takes the posterior distributions from the slices in the peakbag object and samples the joint probability of all the slices.

Notes

This is a coarse estimate of the rotation rate and inclination. A better estimate can be achieved by modeling the entire envelope at once instead of using the slices. To do this set slices=0 when initializing peakbag.

Parameters:

• pb (object) – Peakbagging object containing the instances for analysis.

• NKDE (int, optional) – Number of samples for KDE estimation. Default is 2500.

• bw (float, optional) – Bandwidth for KDE. Default is 0.03.

lnJointPost(theta)

Computes the joint log-posterior probability for the given parameters.

Parameters:

theta (array-like) – Parameter vector.

Returns:

Joint log-posterior probability.

Return type:

lnP float

makeKDEs(N, bw)

Creates KDEs of the posterior distributions from the spectrum slices.

Parameters:

• N (int) – Number of samples for KDE estimation.

• bw (float) – Bandwidth for KDE.

priorLnProb(thetaU)

Computes the log-prior probability for the given parameters.

Parameters:

thetaU (dict) – Dictionary of parameters.

Returns:

lnP – Log-prior probability.

Return type:

float

unpackParams(theta)

Cast the parameters in a dictionary

Parameters:

theta (array) – Array of parameters drawn from the posterior distribution.

Returns:

thetaU – The unpacked parameters.

Return type:

dict

unpackSamples(samples)

Unpacks the sample matrix into a dictionary.

Parameters:

samples (array-like) – Sample matrix.

Returns:

S – Dictionary of samples.

Return type:

dict

class pbjam.peakbagging.peakbag(f, s, ell, freq, height=None, width=None, zeta=None, dnu=None, d02=None, freqLimits=[], rotAsym=None, RV=None, slice=False, dipoleMask=None, snrInput=False, samplerType='emcee', **kwargs)

A class for handling the peakbagging run(s). The spectrum can be divided into a desired number of slices which can be peakbagged independently, which tends to be faster than fitting the whole spectrum at once.

The default setting is to divide the spectrum into roughly the number of radial orders provided in the input frequency list. However, checks are performed to make sure closely spaced modes aren’t split so any potential correlation is correctly accounted for. The number of slices may therefore be less than the requested if it’s not possible to separate the modes.

The slicing is done by a K-means clustering algorithm, where clusters are then merged if they split up for example an l=2,0 mode pair.

Parameters:

• f (array-like) – The frequency array of the spectrum.

• s (array-like) – The values of the power density spectrum.

• ell (array-like) – Angular degree of the modes.

• freq (array-like) – Frequencies of the modes corresponding to the angular degrees.

• height (array-like, optional) – Heights of the modes corresponding to the angular degrees. Default is None, in which case an SNR of 1 is assumed. Providing betters estimates from modeID is however strongly recommended.

• width (array-like, optional) – Widths of the modes corresponding to the angular degrees. Default is None, in which case an width of 0.1 is assumed. Providing betters estimates from modeID is however strongly recommended.

• zeta (array-like, optional) – Mixed-mode coupling factors of the modes corresponding to the angular degrees. Default is None, in which case 0 (no mixing) is assumed for all modes.

• dnu (float, optional) – Large frequency separation. Default is None, in which case it’s estimated from the provided ell and freq arrays.

• d02 (float, optional) – Small frequency separation between l=0 and l=2 modes. Default is None, in which case it’s estimated from the provided ell and freq arrays.

• freqLimits (list, optional) – Frequency limits for the peakbagging. Default is an empty list, in which case the lowest and highest radial mode frequencies -/+ 1 radial order are used to set the limits.

• rotAsym (array-like, optional) – Rotational asymmetry parameters. Default is None.

• RV (array-like, optional) – Radial velocity and associated error of the star in km/s. Default is None.

• slices (int, optional) – Number of slices for mode fitting. Default is -1, in which case the number of radial orders determined from freq is used.

• snrInput (bool, optional) – Flag indicating if the input is signal-to-noise ratio (SNR) spectrum. Default is False.

• **kwargs (dict) – Additional keyword arguments.

N_p

Number of radial (l=0) modes.

Type:

int

Nmodes

Total number of modes.

Type:

int

createPeakbagInstances(sampler)

Create a list of class instances which do the modeling in each spectrum slice. Only the relevant part of the spectrum and mode parameters are passed to each instance.

Unless self.slices is set to 0, attempts to find the best number of slices to use. Starts at the number of l=0 modes, but may adjust down if there are closely spaced sets of modes. Similarly if a choice of self.slices > 0 means a slice has to be placed between two closely spaced modes, the number of slices is also adjusted downward to avoid slicing between closely spaced peaks.

Notes

• If self.slices is 0, no slicing occurs.

• The method stores the created instances in the pbInstances attribute.

getBkg(a=0.66, b=0.88, skips=100)

Estimate the background

Takes an average of the power at linearly spaced points along the log(frequency) axis, where the width of the averaging window increases as a power law.

The mean power values are interpolated back onto the full linear frequency axis to estimate the background noise level at all frequencies.

Returns:

b – Array of psd values approximating the background.

Return type:

array

getRotationInclination()

Computes the joint posterior distribution for rotation and inclination given the samples from the peakbag instances.

Kernel density estimates are computed for all the marginalized posterior distributions. These are then sampled using an MCMC sampler, where the resulting posterior is given by the joint probability of getting the individual posteriors.

Summary statistic and samples are added to the pb.result dictionary.

parseSamples(Nsamples)

Parses the samples to extract and organize the model parameters.

Attempts to include at most N samples from the model, but will default to the actual number of samples of the model parameters if it’s less than N.

The resulting dictionary contains some global parameters, ell, enn, emm etc. and two dictionaries, one containing the samples drawn and one with their summary statistics.

Parameters:

• smp (dict) – A dictionary of sampled parameters.

• Nmax (int, optional) – Maximum number of samples to include. Default is 5000.

Returns:

result – A dictionary containing parsed and organized model parameters.

Return type:

dict

sliceSpectrum()

Slice up the envelope. Sets a series of frequency limits which divide the detected modes into roughly equal number of modes per slice.

Parameters:

• result (dict) – Dictionary of results from the mode ID stage.

• fac (int, optional) – Factor scale dnu to include modes outside envelope, usually very mixed l=1 modes, by default 1.

Returns:

limits – Frequencies delimiting the slices.

Return type:

np.array