Mode Identification¶
PBjam is meant to perform (semi) automated peakbagging of solar-like oscillators. The idea being that it should be relatively easy and painless to get mode frequencies of stars with a precision that is suitable for stellar modeling. Let’s see if we’ve succeeded…
The code can be roughly divided into two main sections, the mode identification stage and the detailed peakbagging stage. This tutorial only covers the mode ID, whereas the peakbagging will be covered in another notebook. The mode ID stage fits a model consisting of a known set of radial orders and angular degrees to the spectrum, so that if a peak is found in the spectrum and the model labels it l=0 for example, then roughly \(\Delta \nu\) further along the spectrum we should expect to see another l=0 mode. If we don’t then we need to vary the parameters of the model, by for example increasing or decreasing \(\Delta \nu\).
This is all done in a Bayesian framework, so we have a prior and we sample the posterior distribution which is the product of the prior and the likelihood of the data given the model. This prior is built from thousands of previous measurements of the model parameters of other stars, so that if we provide an initial estimate of \(\Delta \nu\) and \(\nu_{\mathrm{max}}\) for example, then PBjam knows what range of other mode parameters are important, and crucially also any correlation there might be between them. Of course this all hinges on the solar-like oscillators having a predictable pattern of modes, where it’s just the scale of the pattern that remains the unknown quantity (\(\Delta \nu\), \(\nu_{\mathrm{max}}\) etc.), and so naturally the models in PBjam wouldn’t work on other types of oscillators.
Once a rough pattern is established in the mode ID stage, the mode parameters are passed to the detailed peakbagging stage, where the model is relaxed considerably such that the mode frequencies for example are all left as free parameters. This allows PBjam to capture more detailed structure in the spectrum, which could occur due, for example, to acoustic glitches, or simply to correct for inaccuracies in the mode ID model.
You can read more about the idea behind PBjam in Nielsen et al. (2021) and the current model setup in Nielsen et al. (2025, in prep).
We can start with setting the target name and loading the data. You can use a custom time series or power density spectrum, but the IO class does the tedious work for you, so for now we’ll just use this.
Note: The IO.psd class normalizes the PSD consistently (Parseval). If you are using a custom PSD, be sure it is Parseval normalized, otherwise the built-in prior on mode heights won’t be very helpful. The easiest way to do this is to start with the custom time series and use the IO.psd class to compute the psd, otherwise you will have to normalize the spectrum so the integral of power is equivalent to variance in the light curve.
[10]:
from pbjam import IO
tgt = 'KIC5184732'
psd = IO.psd(tgt, lk_kwargs={'exptime': 60, 'mission':'Kepler', 'author':'Kepler'})
psd()
f = psd.freq
s = psd.powerdensity
For targets with very long time series (~years), and high SNR of the oscillations, it’s usually not necessary to use the spectrum at full resolution. The mode ID process just takes longer. So below we’ve reduced the resolution by a factor of 10, which is still enough to get a good estimate of the background and also the mode ID.
Warning: If you’re working with sub-giants or red giants, some modes can be very narrow and you may therefore need the full resolution.
If you are using a shorter time series, adjust the downsampling factor accordingly. Just be sure you can still make out the modes in the spectrum.
[11]:
downsampling = 10
f = f[::downsampling]
s = s[::downsampling]
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 18
fig, ax = plt.subplots(figsize=(16,9))
ax.plot(f,s)
ax.set_xlim(1500, 2600)
ax.set_ylim(0, 60)
ax.set_ylabel(r'PSD [ppm$^2/\mu$Hz]')
ax.set_xlabel(r'Frequency [$\mu$Hz]')
[11]:
Text(0.5, 0, 'Frequency [$\\mu$Hz]')
Observational parameters¶
As with the previous version of PBjam we supply some basic initial parameters to guide the mode ID process. These are \(T_{\mathrm{eff}}\), Gaia \(G_{\mathrm{bp}}-G_{\mathrm{rp}}\) color index, an estimate of \(\nu_{\mathrm{max}}\) and \(\Delta \nu\). You need to look up these parameters for the target you’ve chosen.
Simbad is very helpful for looking up the \(G_{\mathrm{bp}}-G_{\mathrm{rp}}\) color information, you can use either the GDR2 or GDR3 values here, since they don’t need to be super precise or accurate.
Each bit of information is added to a dictionary, where each entry consists of a list-like entry (tuples or arrays are also fine) which contain the mean and error of the relevant parameter.
[3]:
obs = {'teff' : (5846, 80), # K
'bp_rp': (0.819174, 0.05),
'numax': (2089.3, 20.0), # muHz
'dnu' : (95.545, 0.5)} # muHz
Priors¶
With the update to PBjam we now apply a PCA-based dimensionality reduction method to the sample of prior observations for each of the model parameters. This means that the priors are constructed in a lower-dimensional latent space, instead of the full model space. We then compute 1D KDEs of the latent parameter distributions, which are then used in the sampling.
This relies on the sample of prior observations being fairly dense in the parameter space around the target star. Due to the increased model complexity (\(\ell=1\) modes for example) PBjam doesn’t currently have complete coverage of all parameters in all parts of parameter space. We therefore added an option to include manually specified priors for some of the parameters, in which case they are no longer included in constructing the latent space priors. An example of how to do this is shown at the end of this tutorial. For now let’s get going with the mode ID.
Initializing the mode ID sampler¶
The mode identification in PBjam 2.0 is done in two main steps:
We first identify the \(\ell=2,0\) modes by sampling a simple model of asymptotic p-modes. We also include a background model at this stage. This is almost identical to what was used in the first release of PBjam.
We then divide the PSD by this model to obtain a residual SNR spectrum which mostly just contains \(\ell=1\) modes. We then apply a choice of three different models to this residual to find the \(\ell=1\) modes.
Let’s start by initializing the mode ID and running the \(\ell=2,0\) model. The mode ID class takes the frequency and spectrum arrays we defined above, as well as the dictionary of observational parameters.
As in the first version of PBjam the choice of the number of radial p-mode orders, \(N_p\), is left to the user. Including an excessive amount of orders will probably not improve the result since we rely on the 2nd order asymptotic relation for p-modes. This is only reliable within a few orders of \(\nu_{\mathrm{max}}\) anyway, and so a suitable number for most stars is between 5-10.
[ ]:
from pbjam.modeID import modeID
N_p = 7
M = modeID(f, s, obs, N_p=N_p)
The \(\ell=2,0\) mode ID¶
With the class initialized we can run the \(\ell=2,0\) model. For those who have used PBjam 1.0 this stage should be fairly familiar, however here we are now using Dynesty to do the sampling instead of emcee.
[5]:
M.runl20model();
18863it [04:31, 69.38it/s, +600 | bound: 197 | nc: 1 | ncall: 557270 | eff(%): 3.496 | loglstar: -inf < -7096.259 < inf | logz: -7125.379 +/- 0.196 | dlogz: 0.000 > 0.100]
After running the l20model the model class instance will be updated to contain a samples attribute, which are the likelihood-weighted samples from the nested sampling process. Note that these are samples drawn from the prior that is defined in the latent space, this means that they don’t necessarily correspond to any physical parameters, but are instead linear combinations of these parameters.
Don’t worry though, The mode ID class instance will also be updated to contain the results in a form that is a bit more human readable. These are contained in the result attribute, which is just a big dictionary of summary statistics and samples of the model parameters. The parameters have all been converted to linear scale.
[6]:
M.result
[6]:
{'ell': array([0., 0., 0., 0., 0., 0., 0., 2., 2., 2., 2., 2., 2., 2.]),
'enn': array([17., 18., 19., 20., 21., 22., 23., 16., 17., 18., 19., 20., 21.,
22.]),
'emm': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
'zeta': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
'summary': {'freq': array([[1.75621656e+03, 1.85150156e+03, 1.94694313e+03, 2.04261580e+03,
2.13844582e+03, 2.23440717e+03, 2.33053276e+03, 1.74930398e+03,
1.84459969e+03, 1.94005485e+03, 2.03571457e+03, 2.13155937e+03,
2.22753129e+03, 2.32365879e+03],
[4.55413513e-01, 2.90680532e-01, 1.48138285e-01, 9.99539760e-02,
1.72728390e-01, 2.78185343e-01, 3.69628177e-01, 5.40605848e-01,
3.71327937e-01, 2.20871121e-01, 1.20323316e-01, 1.31073835e-01,
2.07954646e-01, 3.00153762e-01]]),
'height': array([[3.40373183, 5.12217864, 6.36293712, 6.10159078, 4.74017102,
2.99006372, 1.45942334, 2.33236341, 3.55656203, 4.48092988,
4.37493199, 3.44996245, 2.21239206, 1.09984957],
[0.43395854, 0.37763105, 0.41302235, 0.44004288, 0.33133705,
0.1898756 , 0.12195641, 0.31105897, 0.26859343, 0.29480716,
0.31448696, 0.2414944 , 0.13964041, 0.09007912]]),
'width': array([[1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,
1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,
1.47584377, 1.47584377, 1.47584377, 1.47584377],
[0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,
0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,
0.02288559, 0.02288559, 0.02288559, 0.02288559]]),
'rotAsym': array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]),
'dnu': array([95.60011937, 0.16415247]),
'numax': array([1979.76227256, 12.76013212]),
'eps_p': array([1.36608988, 0.0361177 ]),
'd02': array([6.90729371, 0.1032822 ]),
'alpha_p': array([0.00171483, 0.00018255]),
'env_width': array([209.33328481, 9.82147204]),
'env_height': array([3.22559739, 0.21508694]),
'mode_width': array([1.47584377, 0.02288559]),
'teff': array([5928.20304439, 110.14789371]),
'bp_rp': array([0.78161292, 0.01161255]),
'H1_nu': array([1197.57921981, 768.32510693]),
'H1_exp': array([2.04652181, 1.83179639]),
'H_power': array([1397.63266683, 178.4591665 ]),
'H2_nu': array([399.62551993, 52.21881081]),
'H2_exp': array([2.94603882, 0.59286184]),
'H3_power': array([8617.41310727, 4108.05309297]),
'H3_nu': array([0.77586561, 0.45106515]),
'H3_exp': array([4.53051833, 1.3349547 ]),
'shot': array([0.83646182, 0.04979915]),
'nurot_e': array([0.02388869, 0.02063556]),
'inc': array([0.74782077, 0.34498241])},
'samples': {'freq': array([[1755.4437615 , 1851.0341552 , 1946.76619107, ..., 2131.66269706,
2227.81965943, 2324.11826396],
[1756.46359128, 1851.7215938 , 1947.15264206, ..., 2131.72876008,
2227.67894554, 2323.80217672],
[1755.41776398, 1850.99829636, 1946.72201926, ..., 2131.59713308,
2227.75042751, 2324.04691246],
...,
[1756.1572393 , 1851.49169038, 1946.99175917, ..., 2131.50026008,
2227.49718199, 2323.6597216 ],
[1755.3714635 , 1850.95304791, 1946.67543848, ..., 2131.54627802,
2227.69108702, 2323.97670218],
[1756.32998244, 1851.64358947, 1947.12458485, ..., 2131.63882579,
2227.62198624, 2323.77253503]], shape=(5000, 14)),
'height': array([[2.99974512, 4.93146762, 6.32829456, ..., 3.59797248, 2.22161056,
1.06607295],
[3.75432431, 5.35161398, 6.2096966 , ..., 3.27070745, 2.06360308,
1.05510839],
[3.17602407, 5.12527842, 6.50324604, ..., 3.69338995, 2.30383356,
1.12507617],
...,
[3.32420016, 4.9030869 , 5.89878786, ..., 3.3532533 , 2.20826494,
1.18114715],
[3.0440651 , 4.97299756, 6.35429781, ..., 3.6035255 , 2.22872387,
1.0734619 ],
[3.46532403, 5.08369443, 6.09063803, ..., 3.4461285 , 2.26793395,
1.21373934]], shape=(5000, 14)),
'width': array([[1.45672781, 1.45672781, 1.45672781, ..., 1.45672781, 1.45672781,
1.45672781],
[1.48608894, 1.48608894, 1.48608894, ..., 1.48608894, 1.48608894,
1.48608894],
[1.466041 , 1.466041 , 1.466041 , ..., 1.466041 , 1.466041 ,
1.466041 ],
...,
[1.49561911, 1.49561911, 1.49561911, ..., 1.49561911, 1.49561911,
1.49561911],
[1.45491654, 1.45491654, 1.45491654, ..., 1.45491654, 1.45491654,
1.45491654],
[1.49643402, 1.49643402, 1.49643402, ..., 1.49643402, 1.49643402,
1.49643402]], shape=(5000, 14)),
'rotAsym': array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]], shape=(5000, 14)),
'dnu': array([95.87403099, 95.55556281, 95.86492614, ..., 95.64963924,
95.86225478, 95.62965851], shape=(5000,)),
'numax': array([1994.92423044, 1968.12741788, 1993.30710251, ..., 1985.53530824,
1993.94815689, 1984.22900696], shape=(5000,)),
'eps_p': array([1.30527233, 1.37713362, 1.30675012, ..., 1.35531256, 1.30682785,
1.3609634 ], shape=(5000,)),
'd02': array([6.99249224, 6.80511569, 7.00190351, ..., 6.98848979, 6.99636003,
6.94991488], shape=(5000,)),
'alpha_p': array([0.00147738, 0.00181094, 0.00149367, ..., 0.0017315 , 0.00146884,
0.00175038], shape=(5000,)),
'env_width': array([191.98864748, 209.95047023, 194.85381639, ..., 211.12706095,
192.75034356, 211.72628417], shape=(5000,)),
'env_height': array([3.26527321, 3.12038131, 3.34589155, ..., 2.99895515, 3.27415196,
3.09244314], shape=(5000,)),
'mode_width': array([1.45672781, 1.48608894, 1.466041 , ..., 1.49561911, 1.45491654,
1.49643402], shape=(5000,)),
'teff': array([6128.5365841 , 5913.28097489, 6105.14429563, ..., 5958.23734281,
6114.58172715, 5938.6832663 ], shape=(5000,)),
'bp_rp': array([0.76452778, 0.78528912, 0.76774448, ..., 0.77342595, 0.76722184,
0.77885832], shape=(5000,)),
'H1_nu': array([2667.46562177, 1129.54342317, 2627.14469037, ..., 1202.81322221,
2661.52851126, 1208.45154109], shape=(5000,)),
'H1_exp': array([5.60478316, 2.06504413, 5.51922301, ..., 1.96032355, 5.56950369,
1.92396647], shape=(5000,)),
'H_power': array([1102.22232902, 1385.44272929, 1107.69649263, ..., 1415.76550255,
1110.4882708 , 1447.41890278], shape=(5000,)),
'H2_nu': array([497.01580405, 393.09609184, 498.93554013, ..., 398.51846692,
498.77058953, 399.63319227], shape=(5000,)),
'H2_exp': array([1.82594179, 3.0162425 , 1.9919362 , ..., 2.85604736, 1.91246849,
2.8832552 ], shape=(5000,)),
'H3_power': array([ 8627.41302942, 11807.75245977, 7372.96681695, ...,
4767.81585139, 9198.00871568, 5535.88960394], shape=(5000,)),
'H3_nu': array([1.64206018, 0.57638377, 1.38722135, ..., 0.756377 , 1.38039724,
0.72432482], shape=(5000,)),
'H3_exp': array([2.19633475, 4.92352727, 2.06965272, ..., 4.72773144, 2.19553537,
4.61065602], shape=(5000,)),
'shot': array([0.78081294, 0.84982394, 0.76182058, ..., 0.84635884, 0.76220298,
0.85480019], shape=(5000,)),
'nurot_e': array([0.04680386, 0.01347233, 0.05509917, ..., 0.03266987, 0.05369844,
0.02498345], shape=(5000,)),
'inc': array([1.19734749, 0.48226059, 1.23138016, ..., 0.88380619, 1.20072887,
0.65404719], shape=(5000,))}}
The \(\ell=1\) mode ID¶
At this stage, if you only really care about the l=2,0 mode pairs you can skip to the detailed peakbagging. Otherwise we can run the \(\ell=1\) model.
Here the user must make a choice for which model to apply. The choices of models are:
MS: when there are no g-modes near the p-modes. Typically applicable to main-sequence stars. This simply determines the frequencies based on the asymptotic relation for p-modes.SG: when you expect a few g-modes to couple to many p-modes. Typically applicable to sub-giants. Here the mode frequencies are determined by computing the coupling between all combinations of p- and g-modes (see e.g., Ong & Basu 2020).RGB: when you expect lots of g-modes to couple to few p-modes. Typically applicable to low-luminosty red giants. In this model the coupling is computed using the single-value coupling model (e.g. Shibahashi 1979)
Exceptions to the above suggestions of course exist, so if in doubt try different models. All the models should run, but results may vary, and also the time it takes to evaluate the model!
PBjam does have some limited ability to suggest a model, but take care with this as it is just based on drawing some lines in the (\(\Delta\nu\), \(T_{\mathrm{eff}}\)) plane. So when stars start to transition between main-sequence and sub-giants it is bound to be inaccurate.
Note: unless you’re running the MS model, the sampling will probably take a while longer than for the \(\ell=2,0\) since it’s a more complicated model. So consider going for a cup of coffee at this point.
[7]:
M.runl1model();
Input Teff=5846K and dnu=95.545muHz suggests the appropriate l=1 model is: ms
1540it [00:06, 234.58it/s, +150 | bound: 44 | nc: 1 | ncall: 28711 | eff(%): 5.917 | loglstar: -inf < -5958.266 < inf | logz: -5966.166 +/- 0.207 | dlogz: 0.001 > 0.100]
4309it [00:16, 269.31it/s, +150 | bound: 165 | nc: 1 | ncall: 92398 | eff(%): 4.834 | loglstar: -inf < -5958.258 < inf | logz: -5966.166 +/- 0.194 | dlogz: 0.000 > 0.000]
After running the \(\ell=1\) model the result dictionary will be updated to contain the most recent values.
[12]:
M.result
[12]:
{'ell': array([1., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 2., 2., 2.,
2., 2., 2., 2.]),
'enn': array([-1., -1., -1., -1., -1., -1., -1., 17., 18., 19., 20., 21., 22.,
23., 16., 17., 18., 19., 20., 21., 22.]),
'emm': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0.]),
'zeta': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
'summary': {'freq': array([[1.80086678e+03, 1.89615178e+03, 1.99159334e+03, 2.08726602e+03,
2.18309604e+03, 2.27905739e+03, 2.37518298e+03, 1.75621656e+03,
1.85150156e+03, 1.94694313e+03, 2.04261580e+03, 2.13844582e+03,
2.23440717e+03, 2.33053276e+03, 1.74930398e+03, 1.84459969e+03,
1.94005485e+03, 2.03571457e+03, 2.13155937e+03, 2.22753129e+03,
2.32365879e+03],
[9.85093899e-02, 9.85093899e-02, 9.85093899e-02, 9.85093899e-02,
9.85093899e-02, 9.85093899e-02, 9.85093899e-02, 4.55413513e-01,
2.90680532e-01, 1.48138285e-01, 9.99539760e-02, 1.72728390e-01,
2.78185343e-01, 3.69628177e-01, 5.40605848e-01, 3.71327937e-01,
2.20871121e-01, 1.20323316e-01, 1.31073835e-01, 2.07954646e-01,
3.00153762e-01]]),
'height': array([[5.46272976e+00, 7.26705210e+00, 7.85789748e+00, 6.89811267e+00,
4.91045026e+00, 2.83209815e+00, 1.32186648e+00, 3.40373183e+00,
5.12217864e+00, 6.36293712e+00, 6.10159078e+00, 4.74017102e+00,
2.99006372e+00, 1.45942334e+00, 2.33236341e+00, 3.55656203e+00,
4.48092988e+00, 4.37493199e+00, 3.44996245e+00, 2.21239206e+00,
1.09984957e+00],
[2.19692119e-03, 1.36599489e-03, 2.08875927e-04, 1.66699819e-03,
2.24455906e-03, 1.90554424e-03, 1.17507827e-03, 4.33958544e-01,
3.77631051e-01, 4.13022351e-01, 4.40042881e-01, 3.31337051e-01,
1.89875600e-01, 1.21956414e-01, 3.11058973e-01, 2.68593426e-01,
2.94807164e-01, 3.14486955e-01, 2.41494405e-01, 1.39640406e-01,
9.00791207e-02]]),
'width': array([[1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,
1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,
1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,
1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,
1.47584377],
[0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,
0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,
0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,
0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,
0.02288559]]),
'rotAsym': array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0.]]),
'd01': array([44.65021781, 0.09850939]),
'nurot_e': array([0.02388869, 0.02063556]),
'inc': array([0.74782077, 0.34498241]),
'dnu': array([95.60011937, 0.16415247]),
'numax': array([1979.76227256, 12.76013212]),
'eps_p': array([1.36608988, 0.0361177 ]),
'd02': array([6.90729371, 0.1032822 ]),
'alpha_p': array([0.00171483, 0.00018255]),
'env_width': array([209.33328481, 9.82147204]),
'env_height': array([3.22559739, 0.21508694]),
'mode_width': array([1.47584377, 0.02288559]),
'teff': array([5928.20304439, 110.14789371]),
'bp_rp': array([0.78161292, 0.01161255]),
'H1_nu': array([1197.57921981, 768.32510693]),
'H1_exp': array([2.04652181, 1.83179639]),
'H_power': array([1397.63266683, 178.4591665 ]),
'H2_nu': array([399.62551993, 52.21881081]),
'H2_exp': array([2.94603882, 0.59286184]),
'H3_power': array([8617.41310727, 4108.05309297]),
'H3_nu': array([0.77586561, 0.45106515]),
'H3_exp': array([4.53051833, 1.3349547 ]),
'shot': array([0.83646182, 0.04979915])},
'samples': {'freq': array([[1800.91800855, 1896.203005 , 1991.64457353, ..., 2131.66269706,
2227.81965943, 2324.11826396],
[1801.05115263, 1896.33614909, 1991.77771762, ..., 2131.72876008,
2227.67894554, 2323.80217672],
[1800.90252919, 1896.18752564, 1991.62909417, ..., 2131.59713308,
2227.75042751, 2324.04691246],
...,
[1800.90312155, 1896.188118 , 1991.62968653, ..., 2131.46228486,
2227.40020322, 2323.50202836],
[1800.86431731, 1896.14931377, 1991.5908823 , ..., 2131.50285171,
2227.44096778, 2323.55290985],
[1800.74920869, 1896.03420515, 1991.47577368, ..., 2131.46922559,
2227.43713675, 2323.56232498]], shape=(4459, 21)),
'height': array([[5.46387219, 7.26776224, 7.85778856, ..., 3.59797248, 2.22161056,
1.06607295],
[5.46684096, 7.26960619, 7.85750329, ..., 3.27070745, 2.06360308,
1.05510839],
[5.463527 , 7.2675477 , 7.85782153, ..., 3.69338995, 2.30383356,
1.12507617],
...,
[5.46354021, 7.26755591, 7.85782026, ..., 3.12164498, 1.98563627,
1.02281077],
[5.46267484, 7.26701796, 7.85790271, ..., 3.74821676, 2.46381313,
1.32016318],
[5.46010752, 7.26542093, 7.85814568, ..., 3.10471094, 2.01350759,
1.05473943]], shape=(4459, 21)),
'width': array([[1.47584377, 1.47584377, 1.47584377, ..., 1.45672781, 1.45672781,
1.45672781],
[1.47584377, 1.47584377, 1.47584377, ..., 1.48608894, 1.48608894,
1.48608894],
[1.47584377, 1.47584377, 1.47584377, ..., 1.466041 , 1.466041 ,
1.466041 ],
...,
[1.47584377, 1.47584377, 1.47584377, ..., 1.4696376 , 1.4696376 ,
1.4696376 ],
[1.47584377, 1.47584377, 1.47584377, ..., 1.50430534, 1.50430534,
1.50430534],
[1.47584377, 1.47584377, 1.47584377, ..., 1.46569286, 1.46569286,
1.46569286]], shape=(4459, 21)),
'rotAsym': array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]], shape=(4459, 21)),
'd01': array([44.7014464 , 44.83459049, 44.68596704, ..., 44.6865594 ,
44.64775516, 44.53264655], shape=(4459,)),
'nurot_e': array([0.04680386, 0.01347233, 0.05509917, ..., 0.03266987, 0.05369844,
0.02498345], shape=(5000,)),
'inc': array([1.19734749, 0.48226059, 1.23138016, ..., 0.88380619, 1.20072887,
0.65404719], shape=(5000,)),
'dnu': array([95.87403099, 95.55556281, 95.86492614, ..., 95.64963924,
95.86225478, 95.62965851], shape=(5000,)),
'numax': array([1994.92423044, 1968.12741788, 1993.30710251, ..., 1985.53530824,
1993.94815689, 1984.22900696], shape=(5000,)),
'eps_p': array([1.30527233, 1.37713362, 1.30675012, ..., 1.35531256, 1.30682785,
1.3609634 ], shape=(5000,)),
'd02': array([6.99249224, 6.80511569, 7.00190351, ..., 6.98848979, 6.99636003,
6.94991488], shape=(5000,)),
'alpha_p': array([0.00147738, 0.00181094, 0.00149367, ..., 0.0017315 , 0.00146884,
0.00175038], shape=(5000,)),
'env_width': array([191.98864748, 209.95047023, 194.85381639, ..., 211.12706095,
192.75034356, 211.72628417], shape=(5000,)),
'env_height': array([3.26527321, 3.12038131, 3.34589155, ..., 2.99895515, 3.27415196,
3.09244314], shape=(5000,)),
'mode_width': array([1.45672781, 1.48608894, 1.466041 , ..., 1.49561911, 1.45491654,
1.49643402], shape=(5000,)),
'teff': array([6128.5365841 , 5913.28097489, 6105.14429563, ..., 5958.23734281,
6114.58172715, 5938.6832663 ], shape=(5000,)),
'bp_rp': array([0.76452778, 0.78528912, 0.76774448, ..., 0.77342595, 0.76722184,
0.77885832], shape=(5000,)),
'H1_nu': array([2667.46562177, 1129.54342317, 2627.14469037, ..., 1202.81322221,
2661.52851126, 1208.45154109], shape=(5000,)),
'H1_exp': array([5.60478316, 2.06504413, 5.51922301, ..., 1.96032355, 5.56950369,
1.92396647], shape=(5000,)),
'H_power': array([1102.22232902, 1385.44272929, 1107.69649263, ..., 1415.76550255,
1110.4882708 , 1447.41890278], shape=(5000,)),
'H2_nu': array([497.01580405, 393.09609184, 498.93554013, ..., 398.51846692,
498.77058953, 399.63319227], shape=(5000,)),
'H2_exp': array([1.82594179, 3.0162425 , 1.9919362 , ..., 2.85604736, 1.91246849,
2.8832552 ], shape=(5000,)),
'H3_power': array([ 8627.41302942, 11807.75245977, 7372.96681695, ...,
4767.81585139, 9198.00871568, 5535.88960394], shape=(5000,)),
'H3_nu': array([1.64206018, 0.57638377, 1.38722135, ..., 0.756377 , 1.38039724,
0.72432482], shape=(5000,)),
'H3_exp': array([2.19633475, 4.92352727, 2.06965272, ..., 4.72773144, 2.19553537,
4.61065602], shape=(5000,)),
'shot': array([0.78081294, 0.84982394, 0.76182058, ..., 0.84635884, 0.76220298,
0.85480019], shape=(5000,))}}
3.5 Specifying manual prior probability densities in the mode ID stage¶
In the event that your target is in a part of parameter space that is poorly sampled it may sometimes be necessary to provide a prior manually for the mode ID to work.
In PBjam 2.0 this is done by supplying a distribution class instance from the pbjam.distributions module. These function more or less the same as the corresponding scipy.stats versions, with the important distinction that they can be incorporated into a function that is jit-compiled with Jax. Jax is an accelerated linear algebra library that can compile individual functions. Jax is continually being developed, so these functions may become redundant in the future.
An example could be if we are reasonably certain about the p-mode phase offset \(\epsilon_p\), for a star, and so don’t want PBjam to define a prior on its own.
[ ]:
import pbjam.distributions as dist
addPriors = {'eps_p': dist.normal(loc=1.45, scale=0.1)}
M_manual = modeID(f, s, obs, N_p=N_p, addPriors=addPriors)
Before running the sampling it is often useful to do what is known as a prior predictive check. Here samples are drawn from the prior distribution and the model is evaluated. You can then check that the prior you have chosen makes sense. Most of the plotting functions in PBjam will have an option to do this.
[ ]:
M_manual.spectrum(stage='prior')
If this looks vaguely sensible you can try running the sampling for that model.
[ ]:
M_manual.runl20model();
[ ]:
M.echelle()
Warning: You should be careful though since many parameters, like \(\nu_{\mathrm{max}}\) or \(\Delta \nu\), need to be specified in terms of their logarithms. You can check the list of model parameters for the l20 model class by inspecting the variables attribute. Unfortunately this is currently only accessible once the model class has been initialized, this will probably be changed in future updates.