BiCEP Examples Notebook

This notebook runs through examples of how to use BiCEP manually through the BiCEP package. The functionality in this package shows how to get plots out of BiCEP for manual editing with matplotlib. Almost all of the functionality detailed here can be obtained through the GUI, but some users may prefer to write their own code using the standard functions in BiCEP.

Library import

Run the following cell to import the BiCEP module.

import matplotlib.pyplot as plt

Data conversion

We can produce a data file that BiCEP can read internally from the magic format using the generate_arai_plot_table function in the conversion module. We will specify an output name of BiCEP_converted.csv. This gets dumped in the folder with the datafiles. For your convenience, we have included a set of MagIC formatted data to convert to the BiCEP format. This may take a few minutes

from BiCEP.conversion import generate_arai_plot_table
project='example_magic_files/'
fname = 'BiCEP_converted'
generate_arai_plot_table(fname,project)

Data conversion finished- check thellier_convert.log for errors

When imported, BiCEP data is stored internally as a set of nested objects, one for the dataset which contains a object for each site (a set of samples expected to have been magnetized in the same field), or each sample (an object collected in the field) which are described using a SpecimenCollection object. Each of these SpecimenCollection objects contains an object for each specimen, which is a subsample that a set of independent measurements are made on. To make the data class, we can initialize a ThellierData object using our BiCEP_converted.csv file. This may take a few seconds to create all of the objects. Printing our ThellierData object allows us to see the number of sites and specimens.

from BiCEP.datamodel import ThellierData
data=ThellierData(project+fname+'.csv')
print(data)

Set of Thellier Data Containing the Sites:
1991-1992 Eruption Site (53 specimens)
hw108   (23 specimens)
hw123   (12 specimens)
hw126   (13 specimens)
hw128   (26 specimens)
hw201   (12 specimens)
hw226   (11 specimens)
hw241   (18 specimens)
BR06    (3 specimens)
P       (36 specimens)
VM      (18 specimens)
BBQ     (12 specimens)
rs25    (5 specimens)
rs26    (5 specimens)
rs27    (10 specimens)
remag-rs61      (6 specimens)
remag-rs62      (6 specimens)
remag-rs63      (5 specimens)
remag-rs78      (4 specimens)
kf      (3 specimens)
Hawaii 1960 Flow        (22 specimens)
SW      (19 specimens)
TS      (53 specimens)
ET1     (3 specimens)
ET2     (2 specimens)
ET3     (4 specimens)
Synthetic60     (7 specimens)
LV      (45 specimens)
MSH     (19 specimens)
FreshTRM        (24 specimens)

Data processing and analysis

Looking at site and specimen data

We can index for our sites and specimens using a key for the site name. Let’s look at site “hw126” and print out the object to see the specimens. This site is a lava flow erupted in Hawaii in 1935 and originally analyzed by Cromwell et al (2015). Because this eruption occurred in the age of instrumentation, we know that the original field this lava flow was erupted in was approximately 36.4 μT.

site=data['hw126']
print(site)

Site containing the specimens:
hw126a1
hw126a2
hw126a3
hw126a4
hw126a5
hw126b1
hw126b2
hw126b3
hw126b4
hw126b5
hw126a6
hw126a7
hw126a8

Let’s take a look at some of these data from one of the specimens from this site. We can use the plot_arai() method on our specimen to plot the Thellier data on an Arai plot (Nagata & Arai, 1963), which shows the magnitude of the magnetization lost by cooling in zero field against the magnetization gained when cooling in a known field. We can also plot a Zijderveld plot using the plot_zijd method, which plots the magnetization vectors of the zero field steps. In the zijderveld plot, black circles represent the component of the magnetization in the x,y plane, whereas red squares represent the component in the x,z plane.

specimen=site['hw126a1']
fig,ax=plt.subplots(1,2,figsize=(12,4))
specimen.plot_arai(ax[0])
specimen.plot_zijd(ax[1])
ax[0].set_title('Arai Plot')
ax[1].set_title('Zijderveld Plot');

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

Both our Arai plot and Zijderveld plot are a straight line. Because we’re plotting magnetization lost in zero field, against magnetization gained in a known field, the ratio of our two lines is the ratio of the Ancient field to the Lab field which is a little under 2. We can check the lab field using the B_lab attribute of the specimen.

specimen.B_lab

np.float64(20.0)

The lab field was 20 μT. If we expect a value of 36.4 μT, then a slope of about 1.8 is what we might expect. Let’s look at another specimen.

specimen=site['hw126a2']
fig,ax=plt.subplots(1,2,figsize=(12,4))
specimen.plot_arai(ax[0])
specimen.plot_zijd(ax[1])
ax[0].set_title('Arai Plot')
ax[1].set_title('Zijderveld Plot');

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

The data for this specimen do not follow a straight line. Although the line slope appears to start out at a reasonable value, it quickly becomes shallower at higher temperatures. If we were to fit a line to these data we might end up with an underestimate of the paleointensity, and there is no reason to suggest which temperature range to pick for our intensity data. Some strict selection criteria might exclude this specimen, others might include it, but could include an incorrect temperature range, which could give a biased paleointensity estimate.

The BiCEP method- accounting for bias

BiCEP works differently. It fits a circle to a scaled set of paleointensity data and uses the tangent joining the circle center and the origin. BiCEP generates many circle fits from a probability distribution. Let’s look at the circle fits that BiCEP samples from the posterior distribution, and the median tangent which we use to fit the line. BiCEP fits are performed at a site level for a reason which will become apparent later. We can use the method site.BiCEP_fit() to fit the BiCEP method (this may take several minutes depending on the speed of your processor and memory available) and plot the circle and tangent fits to the data.

site.BiCEP_fit()
#This cell may generate warnings from the internal pystan code.
#The main error to watch for is one about R_hat, which indicates the sampler hasn't converged, a "bad" run

19:21:10 - cmdstanpy - INFO - CmdStan start processing

19:21:44 - cmdstanpy - INFO - CmdStan done processing.

19:21:45 - cmdstanpy - WARNING - Some chains may have failed to converge.
        Chain 1 had 138 divergent transitions (0.5%)
        Chain 2 had 92 divergent transitions (0.3%)
        Chain 3 had 178 divergent transitions (0.6%)
        Chain 4 had 171 divergent transitions (0.6%)
        Use the "diagnose()" method on the CmdStanMCMC object to see further information.

We can use the specimen.plot_circ() method to draw a circle under the Arai plot data. The tangent=True argument shows the tangent to the circle used to calculate the paleointensity for this specimen.

fig,ax=plt.subplots(1,2,figsize=(12,4))
specimen.plot_circ(ax[0],tangent=True)
specimen.plot_arai(ax[0])
specimen.plot_zijd(ax[1])
ax[0].set_title('Arai Plot')
ax[1].set_title('Zijderveld Plot');

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

We notice that this tangent is slightly less than the initial slope, how does this help us?

BiCEP assumes that the curvature, \(\vec{k}\) (1/radius with sign depending on direction of curvature, Paterson, 2011) of the circle fit is proportional to the calculated ancient field (\(B_{\textrm{anc}}\)). For \(\vec{k}>0\), our “circle” has infinite radius, and so becomes a perfect line. BiCEP tries to fit a linear relationship between \(B_{\textrm{anc}}\) and \(\vec{k}\) for all specimens in a site. The value can be extrapolated (or interpolated) to \(\vec{k}\) = 0 to give us our unbiased paleointensity estimate. In the case of the plotted specimen hw126a2, BiCEP predicts that this specimen should underestimate the site paleointensity slightly, due to its curvature.

We can plot the relationship between and using the method site.regplot(). We plot the expected field value as a red line.

fig,ax=plt.subplots()
site.regplot(ax)

#x limits, y limits and y label are not set by this method as it is used slightly differently in our GUI.

ax.set_xlim(-0.2,1.2)
ax.set_ylim(15,50);
ax.set_ylabel('$B_{anc}$ ($\mu$T)')
ax.axhline(36.4,color='r',zorder=0);

The green dots with black errorbars are our individual specimen estimates for \(\vec{k}\) and \(B_{\textrm{anc}}\) and their 95% confidence intervals. The blue lines are our BiCEP fits to these data. Notice that these lines cross the \(\vec{k}=0\) line close to the expected field value, yielding an accurate and precise estimate of the paleointensity, without excluding any specimens! We can get a sense of the distribution of estimates for using the site.histplot() method.

fig,ax=plt.subplots()
site.histplot(ax)

ax.set_xlim(30,40)
ax.axvline(36.4,color='r'); #Actual field at this site

We can see that our estimate for \(B_{\textrm{anc}}\) is close to the median value of this distribution, and well within the 95% confidence interval denoted by the black line. Our estimate is accurate and precise, with a 95% confidence interval of ~4 μT.

Data manipulation and choosing interpretations.

You might notice that some of our specimens have larger error bars on their \(\vec{k}\) estimates than others in our BiCEP fitting figure. Let’s investigate what’s going on with these specimens.

fig,ax=plt.subplots(1,2,figsize=(12,4));
hw126a7=site['hw126a7']
hw126a8=site['hw126a8']
hw126a7.plot_circ(ax[0],tangent=True)
hw126a7.plot_arai(ax[0])
ax[0].set_title('hw126a7')
hw126a8.plot_circ(ax[1],tangent=True)
hw126a8.plot_arai(ax[1])
ax[1].set_title('hw126a8');

For these specimens, our analysis might be slightly different. Our specimen “hw126a7” has a pTRM check (triangle) which is very displaced from the line at the 300 °C temperature step. This is a repeat of the 300 °C in-field measurement after heating to 350 °C in zero field. It indicates a chemical alteration of the specimen between 300 °C and 350 °C. We can use this by calculating the DRAT statistic (Selkin & Tauxe, 2000), using the attribute specimen.drat, and noticing that it is high (criteria vary, but 15 is considered a high value by any metric). Although we have already demonstrated the effectiveness of BiCEP without using this statistic, we have evidence that we should only be using the temperatures up to 300 °C for this specimen. We can do this using the method specimen.change_temps(), and we save this change using the method specimen.save_changes(). Note that our attributes, such as specimen.drat are affected by the change in temperatures.

Our specimen “hw126a8” has a very early failed pTRM check, and fitting a circle seems like it would probably be inappropriate. We can elect to exclude this specimen as no part of the line is likely to give a good paleointensity estimate. This can be done by setting the attribute specimen.active=False Ultimately, due to its large uncertainty in the circle fit, it has very little effect on the estimate of \(B_{\textrm{anc}}\) in the first place.

The MAD statistic of Kirshvink (1980) and the DANG Statistic (Tanaka & Kobayashi, 2003, Tauxe & Staudigel, 2004) are also available using the methods specimen.mad and specimen.dang.

#Printing the DRAT parameter of the specimen to
print('hw126a7 DRAT: %2.1f'%hw126a7.drat)
print('hw126a8 DRAT: %2.1f'%hw126a8.drat)

lowerTemp=0
upperTemp=300
lowerTempK,upperTempK=lowerTemp+273,upperTemp+273

hw126a7.change_temps(lowerTempK,upperTempK)
hw126a7.save_changes()
hw126a8.active=False

hw126a7 DRAT: 0.0
hw126a8 DRAT: 0.0

We now run the BiCEP method again with these changes in place, and plot up our site results with the results for these two specimens.

site.BiCEP_fit()

19:21:47 - cmdstanpy - INFO - CmdStan start processing

19:22:23 - cmdstanpy - INFO - CmdStan done processing.

19:22:23 - cmdstanpy - WARNING - Some chains may have failed to converge.
        Chain 1 had 684 divergent transitions (2.3%)
        Chain 2 had 213 divergent transitions (0.7%)
        Chain 3 had 889 divergent transitions (3.0%)
        Chain 4 had 838 divergent transitions (2.8%)
        Use the "diagnose()" method on the CmdStanMCMC object to see further information.

#set up subplots
fig,ax=plt.subplots(2,2,figsize=(12,8))

#Plot our Arai plots with circle fits
hw126a7.plot_circ(ax[0,0])
hw126a7.plot_arai(ax[0,0])

hw126a8.plot_circ(ax[0,1])
hw126a8.plot_arai(ax[0,1])


#Plot our plot of the BiCEP fit.
site.regplot(ax[1,0])
ax[1,0].set_xlim(-0.2,1.2)
ax[1,0].set_ylim(15,50);
ax[1,0].set_ylabel('$B_{anc}$ ($\mu$T)')
ax[1,0].axhline(36.4,color='r',zorder=0);


#Plot our histogram
site.histplot(ax[1,1])
ax[1,1].set_xlim(30,40)
ax[1,1].axvline(36.4,color='r');

We notice from our histogram that doing this data manipulation has very slightly improved our accuracy and precision- but the overall distribution has not changed vary much. We notice from our Arai plot for specimen hw126a7 that the uncertainty in our circle fits has become much larger by excluding these temperature steps. This is because we assume the uncertainty in the circle fit for the whole Arai plot, rather than just the segment we’re looking at. In this way, we are penalized for excluding data from the Arai plot, without having to use some kind of “length of the line” criterion.

Auto Interpreting With BiCEP

Specimens can be auto-interpreted using the auto_interpret function from the criteria module. This function allows us to auto-interpret by first applying minimal criteria, (MAD, DANG and DRAT) to ensure that we have an interpretation which is origin-trending and for which there is no alteration. Within these specimens, we then take the interpretation that maximizes FRAC.

We only use specimens for which DANG and DRAT are <10 and MAD \(_{\text{Coe}}\)<5 and a minimum number of steps of four. For more information see Paterson et al (2014) . For information on MAD \(_{\text{Coe}}\), please refer to Cych et. al. (2023).

The following functions show how one can auto interpret data using BiCEP. You do not need to rerun these analyses to plot the figures as the BiCEP fits are saved in .nc files in this folder, and the interpretation temperatures are saved to a .redo file.

from BiCEP.criteria import auto_interpret
#Example of auto_interpreting BiCEP_data
#Get samples where MAD_coe < 5, DANG < 10, DRAT < 10
auto_interpret(site,mad=5,dang=10,drat=10,mad_type='mad_coe')
site.BiCEP_fit()

19:22:30 - cmdstanpy - INFO - CmdStan start processing

19:23:21 - cmdstanpy - INFO - CmdStan done processing.

19:23:21 - cmdstanpy - WARNING - Some chains may have failed to converge.
        Chain 1 had 614 divergent transitions (2.0%)
        Chain 2 had 953 divergent transitions (3.2%)
        Chain 3 had 838 divergent transitions (2.8%)
        Chain 4 had 1639 divergent transitions (5.5%)
        Use the "diagnose()" method on the CmdStanMCMC object to see further information.

#set up subplots
fig,ax=plt.subplots(2,2,figsize=(12,8))

#Plot our Arai plots with circle fits
hw126a7.plot_circ(ax[0,0])
hw126a7.plot_arai(ax[0,0])

hw126a8.plot_circ(ax[0,1])
hw126a8.plot_arai(ax[0,1])


#Plot our plot of the BiCEP fit.
site.regplot(ax[1,0])
ax[1,0].set_xlim(-0.2,1.2)
ax[1,0].set_ylim(15,50);
ax[1,0].set_ylabel('$B_{anc}$ ($\mu$T)')
ax[1,0].axhline(36.4,color='r',zorder=0);


#Plot our histogram
site.histplot(ax[1,1])
ax[1,1].set_xlim(30,40)
ax[1,1].axvline(36.4,color='r');

Notice that to satisfy MAD and DRAT, even fewer points from specimen hw126a7 were selected. This has led to the enormous error bars on the estimate from this specimen, seen in the bottom left figure. Sometimes this causes problems with the model fit if the number of points is too small. If this is the case, then try setting the specimen.active attribute to False for that specimen.

Automating, diagnostics and saving fits

We can test whether a fit has yielded a successful paleointensity using the get_sampler_succcess() function from the extract module. Check the docstring to see options - but by default this only accepts things which are an A or a B grade.

from BiCEP.extract import get_sampler_success

get_sampler_success(site)

hw126 successful! Estimated intensity 35.8 with bounds of 34.1 - 37.3 μT

True

We can automate the entire process of getting obtaining a paleointensity and checking whether it’s successful using the run_site_fit() function from the extract module. The save_nc argument can be used to save to netcdf if the result passes.

from BiCEP.extract import run_site_fit
run_site_fit(site,save_nc=True)

19:23:39 - cmdstanpy - INFO - CmdStan start processing

19:23:49 - cmdstanpy - INFO - CmdStan done processing.

19:23:50 - cmdstanpy - WARNING - Some chains may have failed to converge.
        Chain 1 had 79 divergent transitions (0.8%)
        Chain 2 had 55 divergent transitions (0.5%)
        Chain 3 had 21 divergent transitions (0.2%)
        Chain 4 had 69 divergent transitions (0.7%)
        Use the "diagnose()" method on the CmdStanMCMC object to see further information.

hw126 successful! Estimated intensity 35.9 with bounds of 34.3 - 37.4 μT

You should see a file called hw126.nc has appeared in the folder. This is our saved fit that can be reopened by the GUI next time.

Advanced user options - tweaks and accessing raw data

The “fit” object

For each site, BiCEP creates a fit attribute which contains all the information for the fit to that specimen. This is an arviz.InferenceData object

site.fit

arviz.InferenceData

• posterior

  <xarray.Dataset> Size: 23MB
  Dimensions:             (chain: 4, draw: 10000, k_dim_0: 10,
                           dist_to_edge_dim_0: 10, int_real_dim_0: 10,
                           phi_dim_0: 10, R_dim_0: 10, x_c_dim_0: 10,
                           y_c_dim_0: 10)
  Coordinates:
    * chain               (chain) int64 32B 0 1 2 3
    * draw                (draw) int64 80kB 0 1 2 3 4 ... 9995 9996 9997 9998 9999
    * k_dim_0             (k_dim_0) int64 80B 0 1 2 3 4 5 6 7 8 9
    * dist_to_edge_dim_0  (dist_to_edge_dim_0) int64 80B 0 1 2 3 4 5 6 7 8 9
    * int_real_dim_0      (int_real_dim_0) int64 80B 0 1 2 3 4 5 6 7 8 9
    * phi_dim_0           (phi_dim_0) int64 80B 0 1 2 3 4 5 6 7 8 9
    * R_dim_0             (R_dim_0) int64 80B 0 1 2 3 4 5 6 7 8 9
    * x_c_dim_0           (x_c_dim_0) int64 80B 0 1 2 3 4 5 6 7 8 9
    * y_c_dim_0           (y_c_dim_0) int64 80B 0 1 2 3 4 5 6 7 8 9
  Data variables:
      int_site            (chain, draw) float64 320kB 36.5 34.2 ... 36.33 37.26
      sd_site             (chain, draw) float64 320kB 0.7572 1.624 ... 1.103 2.213
      k                   (chain, draw, k_dim_0) float64 3MB 0.01835 ... -0.1276
      dist_to_edge        (chain, draw, dist_to_edge_dim_0) float64 3MB 0.6975 ...
      c                   (chain, draw) float64 320kB -16.08 -18.26 ... -23.17
      int_real            (chain, draw, int_real_dim_0) float64 3MB 36.86 ... 3...
      phi                 (chain, draw, phi_dim_0) float64 3MB 0.7863 ... 0.796
      R                   (chain, draw, R_dim_0) float64 3MB 54.5 2.967 ... -7.838
      x_c                 (chain, draw, x_c_dim_0) float64 3MB 39.0 ... -4.96
      y_c                 (chain, draw, y_c_dim_0) float64 3MB 39.07 ... -5.066
  Attributes:
      created_at:                 2025-08-12T18:23:50.607928+00:00
      arviz_version:              0.21.0
      inference_library:          cmdstanpy
      inference_library_version:  1.2.5

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 10000
    • k_dim_0: 10
    • dist_to_edge_dim_0: 10
    • int_real_dim_0: 10
    • phi_dim_0: 10
    • R_dim_0: 10
    • x_c_dim_0: 10
    • y_c_dim_0: 10
  • Coordinates: (9)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 9996 9997 9998 9999

      array([   0,    1,    2, ..., 9997, 9998, 9999], shape=(10000,))

    • k_dim_0
      (k_dim_0)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

    • dist_to_edge_dim_0
      (dist_to_edge_dim_0)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

    • int_real_dim_0
      (int_real_dim_0)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

    • phi_dim_0
      (phi_dim_0)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

    • R_dim_0
      (R_dim_0)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

    • x_c_dim_0
      (x_c_dim_0)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

    • y_c_dim_0
      (y_c_dim_0)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

  • Data variables: (10)
    • int_site
      (chain, draw)
      float64
      36.5 34.2 36.32 ... 36.33 37.26

      array([[36.4977, 34.2017, 36.3226, ..., 35.7798, 36.0096, 35.4299],
             [36.827 , 36.4198, 34.0576, ..., 35.3379, 35.7266, 36.4886],
             [35.493 , 36.6538, 35.79  , ..., 35.6682, 35.9937, 36.4254],
             [36.6561, 35.4657, 37.0806, ..., 35.1644, 36.3285, 37.2612]],
            shape=(4, 10000))

    • sd_site
      (chain, draw)
      float64
      0.7572 1.624 3.225 ... 1.103 2.213

      array([[0.757206, 1.62438 , 3.22452 , ..., 1.36494 , 1.90146 , 1.86815 ],
             [1.45887 , 1.41425 , 1.52028 , ..., 1.74294 , 2.35597 , 1.82271 ],
             [2.23533 , 0.900743, 1.25724 , ..., 1.16217 , 1.7356  , 0.553311],
             [3.43942 , 2.42982 , 1.48544 , ..., 1.55284 , 1.10321 , 2.21313 ]],
            shape=(4, 10000))

    • k
      (chain, draw, k_dim_0)
      float64
      0.01835 0.337 ... -0.06641 -0.1276

      array([[[ 1.83477e-02,  3.37019e-01,  4.32233e-01, ...,  7.17499e-03,
                5.46733e-04, -1.12214e-01],
              [ 7.72085e-02,  3.18349e-01,  3.81747e-01, ..., -3.18338e-02,
               -1.20265e-01, -1.45062e-01],
              [ 3.68983e-02,  3.76259e-01,  2.76920e-01, ..., -7.27898e-03,
               -1.92277e-02, -1.56140e-01],
              ...,
              [ 1.95842e-02,  2.82486e-01,  3.57004e-01, ..., -6.29515e-02,
               -1.43722e-01, -1.65220e-01],
              [-1.79720e-02,  3.78225e-01,  2.24665e-01, ..., -1.89587e-02,
               -1.28820e-01, -1.36761e-01],
              [ 3.58457e-02,  3.74074e-01,  3.83742e-01, ...,  7.92929e-03,
               -1.20870e-01, -2.11174e-01]],
      
             [[ 3.75149e-02,  3.05107e-01,  3.17356e-01, ..., -1.80268e-02,
               -5.63212e-02, -2.91254e-02],
              [ 4.08837e-03,  3.26285e-01,  3.79208e-01, ..., -6.17128e-03,
               -6.18114e-02, -1.05549e-01],
              [ 3.23776e-02,  3.17421e-01,  3.91406e-01, ..., -2.76009e-02,
               -1.15427e-01, -1.93895e-01],
      ...
              [-9.54077e-03,  3.18113e-01,  4.21588e-01, ..., -3.88341e-02,
               -1.11022e-01, -1.44113e-01],
              [ 5.47108e-02,  3.47408e-01,  3.56644e-01, ..., -1.75755e-02,
               -6.55450e-02, -1.04757e-01],
              [-2.75972e-02,  3.47177e-01,  4.64920e-01, ..., -5.51888e-03,
               -4.94127e-02, -5.61858e-02]],
      
             [[ 5.28950e-02,  3.43898e-01,  3.85070e-01, ...,  6.29879e-04,
               -4.52535e-02, -1.40333e-01],
              [ 1.31691e-01,  3.48224e-01,  3.37088e-01, ..., -2.57052e-02,
               -6.86586e-02, -1.69372e-01],
              [ 1.56000e-02,  3.63638e-01,  4.04150e-01, ...,  1.22358e-02,
               -9.07226e-02, -2.52137e-02],
              ...,
              [ 4.32213e-02,  3.28507e-01,  2.76834e-01, ...,  3.73611e-02,
               -9.21663e-02, -1.41003e-01],
              [ 1.85485e-02,  3.14668e-01,  5.58641e-01, ..., -5.51346e-02,
               -6.02712e-02, -7.53856e-02],
              [ 3.73359e-02,  3.70368e-01,  2.98420e-01, ...,  7.33024e-02,
               -6.64097e-02, -1.27585e-01]]], shape=(4, 10000, 10))

    • dist_to_edge
      (chain, draw, dist_to_edge_dim_0)
      float64
      0.6975 0.6035 ... 0.7202 0.7482

      array([[[0.697545, 0.603518, 0.590005, ..., 0.703924, 0.705586,
               0.741424],
              [0.679517, 0.604079, 0.591771, ..., 0.708039, 0.734576,
               0.738287],
              [0.689874, 0.602198, 0.607157, ..., 0.711001, 0.706012,
               0.744426],
              ...,
              [0.689112, 0.610823, 0.604755, ..., 0.718039, 0.736715,
               0.746308],
              [0.697521, 0.600785, 0.613927, ..., 0.710984, 0.745861,
               0.742974],
              [0.684814, 0.598865, 0.595792, ..., 0.701235, 0.730013,
               0.761007]],
      
             [[0.68386 , 0.604364, 0.592533, ..., 0.712623, 0.722919,
               0.732228],
              [0.688173, 0.606436, 0.598703, ..., 0.708705, 0.721467,
               0.747256],
              [0.683451, 0.603843, 0.589022, ..., 0.707826, 0.72303 ,
               0.743746],
      ...
              [0.698115, 0.601414, 0.593366, ..., 0.710051, 0.725787,
               0.730221],
              [0.687062, 0.608421, 0.596282, ..., 0.703528, 0.720476,
               0.745906],
              [0.702626, 0.605697, 0.592113, ..., 0.713964, 0.713278,
               0.72822 ]],
      
             [[0.692174, 0.605549, 0.580272, ..., 0.707702, 0.71591 ,
               0.741368],
              [0.673901, 0.598074, 0.604281, ..., 0.714309, 0.722642,
               0.75243 ],
              [0.700453, 0.602409, 0.579061, ..., 0.710519, 0.718565,
               0.715625],
              ...,
              [0.686993, 0.60003 , 0.603324, ..., 0.700508, 0.723172,
               0.75143 ],
              [0.692581, 0.607661, 0.580684, ..., 0.713998, 0.716477,
               0.726811],
              [0.686158, 0.60384 , 0.608105, ..., 0.718525, 0.720153,
               0.748222]]], shape=(4, 10000, 10))

    • c
      (chain, draw)
      float64
      -16.08 -18.26 ... -15.1 -23.17

      array([[-16.0812 , -18.2578 ,  -2.6955 , ..., -17.3378 , -11.2183 ,
              -15.9182 ],
             [-17.8348 , -24.3989 , -13.0336 , ..., -11.6405 , -13.2784 ,
              -21.2234 ],
             [-13.5699 , -15.8908 , -15.7251 , ..., -15.0354 , -16.2262 ,
              -15.7836 ],
             [-22.5969 ,  -9.78089, -20.358  , ..., -19.3579 , -15.1032 ,
              -23.1681 ]], shape=(4, 10000))

    • int_real
      (chain, draw, int_real_dim_0)
      float64
      36.86 32.1 28.4 ... 37.36 38.29

      array([[[36.8602, 32.0997, 28.402 , ..., 36.1545, 38.3043, 38.1164],
              [36.3435, 31.5962, 27.9871, ..., 36.0161, 37.663 , 36.298 ],
              [37.1541, 32.9084, 27.529 , ..., 36.1242, 37.2094, 37.3846],
              ...,
              [37.6268, 31.9365, 30.1975, ..., 36.4193, 37.9501, 38.4136],
              [36.1763, 32.4391, 27.7848, ..., 36.3325, 38.2308, 38.8194],
              [35.5684, 31.8202, 27.6918, ..., 35.7576, 38.6205, 38.6714]],
      
             [[36.0984, 31.9778, 27.3209, ..., 36.5403, 37.5817, 38.2747],
              [35.602 , 32.0422, 27.8606, ..., 36.8534, 37.8977, 40.2096],
              [36.4715, 31.9208, 27.548 , ..., 34.7791, 36.8485, 35.3314],
              ...,
              [36.3855, 31.8415, 28.2424, ..., 36.0315, 38.192 , 38.4678],
              [36.2597, 32.427 , 26.9053, ..., 36.2087, 37.4771, 38.3734],
              [37.3029, 32.0992, 28.0352, ..., 36.2394, 37.8821, 37.4181]],
      
             [[35.9644, 32.2836, 27.642 , ..., 36.6559, 36.5029, 37.057 ],
              [37.3574, 32.0924, 27.7501, ..., 36.1247, 38.5478, 38.1156],
              [35.6661, 32.2955, 27.3873, ..., 36.472 , 36.6574, 36.3941],
              ...,
              [36.4275, 32.0273, 28.1212, ..., 36.3772, 37.3182, 35.9966],
              [36.5407, 32.4078, 27.39  , ..., 36.1397, 37.8817, 39.5556],
              [36.4633, 31.6944, 28.9387, ..., 36.3524, 37.2446, 37.751 ]],
      
             [[37.1221, 32.128 , 27.5198, ..., 36.1778, 36.8092, 37.627 ],
              [36.0156, 32.2748, 28.088 , ..., 36.1854, 38.7618, 37.736 ],
              [36.6913, 31.959 , 27.281 , ..., 36.3761, 36.6721, 37.7472],
              ...,
              [36.5329, 31.6408, 27.8299, ..., 36.1132, 37.35  , 38.2398],
              [36.8949, 32.1588, 28.1992, ..., 36.1054, 36.9569, 37.0871],
              [36.5378, 32.1608, 28.1264, ..., 36.7244, 37.3606, 38.2931]]],
            shape=(4, 10000, 10))

    • phi
      (chain, draw, phi_dim_0)
      float64
      0.7863 0.7721 ... 0.8048 0.796

      array([[[0.786275, 0.772119, 0.727709, ..., 0.789153, 0.792326,
               0.798335],
              [0.793333, 0.780023, 0.735024, ..., 0.79107 , 0.800767,
               0.822743],
              [0.782304, 0.759688, 0.743241, ..., 0.789572, 0.806822,
               0.808021],
              ...,
              [0.775984, 0.774667, 0.69739 , ..., 0.785504, 0.796971,
               0.794452],
              [0.795639, 0.766863, 0.738635, ..., 0.786697, 0.793288,
               0.789199],
              [0.804108, 0.77649 , 0.740304, ..., 0.794672, 0.788216,
               0.791108]],
      
             [[0.796715, 0.774022, 0.747023, ..., 0.783846, 0.801847,
               0.796263],
              [0.803637, 0.773015, 0.737279, ..., 0.77958 , 0.797662,
               0.771608],
              [0.791575, 0.774912, 0.742898, ..., 0.808538, 0.811688,
               0.836186],
      ...
              [0.79218 , 0.773247, 0.732647, ..., 0.786083, 0.805363,
               0.826898],
              [0.790628, 0.767345, 0.745763, ..., 0.789358, 0.797872,
               0.779806],
              [0.791689, 0.77847 , 0.718422, ..., 0.786424, 0.806348,
               0.803149]],
      
             [[0.782735, 0.771679, 0.743408, ..., 0.788832, 0.812222,
               0.804792],
              [0.797864, 0.7694  , 0.733235, ..., 0.788726, 0.786391,
               0.803348],
              [0.788571, 0.774315, 0.747751, ..., 0.786098, 0.814084,
               0.8032  ],
              ...,
              [0.790735, 0.779317, 0.737828, ..., 0.789724, 0.804936,
               0.79672 ],
              [0.785804, 0.7712  , 0.73127 , ..., 0.789833, 0.810222,
               0.812012],
              [0.790668, 0.771168, 0.732555, ..., 0.781334, 0.804796,
               0.796024]]], shape=(4, 10000, 10))

    • R
      (chain, draw, R_dim_0)
      float64
      54.5 2.967 2.314 ... -15.06 -7.838

      array([[[ 5.45026e+01,  2.96719e+00,  2.31357e+00, ...,  1.39373e+02,
                1.82905e+03, -8.91152e+00],
              [ 1.29519e+01,  3.14121e+00,  2.61954e+00, ..., -3.14132e+01,
               -8.31500e+00, -6.89359e+00],
              [ 2.71016e+01,  2.65774e+00,  3.61116e+00, ..., -1.37382e+02,
               -5.20084e+01, -6.40452e+00],
              ...,
              [ 5.10616e+01,  3.54000e+00,  2.80109e+00, ..., -1.58852e+01,
               -6.95789e+00, -6.05253e+00],
              [-5.56422e+01,  2.64393e+00,  4.45107e+00, ..., -5.27463e+01,
               -7.76277e+00, -7.31202e+00],
              [ 2.78973e+01,  2.67327e+00,  2.60592e+00, ...,  1.26115e+02,
               -8.27337e+00, -4.73543e+00]],
      
             [[ 2.66561e+01,  3.27754e+00,  3.15104e+00, ..., -5.54728e+01,
               -1.77553e+01, -3.43343e+01],
              [ 2.44597e+02,  3.06480e+00,  2.63707e+00, ..., -1.62041e+02,
               -1.61782e+01, -9.47431e+00],
              [ 3.08855e+01,  3.15039e+00,  2.55489e+00, ..., -3.62307e+01,
               -8.66351e+00, -5.15743e+00],
      ...
              [-1.04813e+02,  3.14354e+00,  2.37199e+00, ..., -2.57506e+01,
               -9.00719e+00, -6.93901e+00],
              [ 1.82779e+01,  2.87846e+00,  2.80391e+00, ..., -5.68972e+01,
               -1.52567e+01, -9.54587e+00],
              [-3.62356e+01,  2.88037e+00,  2.15091e+00, ..., -1.81196e+02,
               -2.02377e+01, -1.77981e+01]],
      
             [[ 1.89054e+01,  2.90784e+00,  2.59693e+00, ...,  1.58761e+03,
               -2.20977e+01, -7.12592e+00],
              [ 7.59351e+00,  2.87171e+00,  2.96658e+00, ..., -3.89026e+01,
               -1.45648e+01, -5.90415e+00],
              [ 6.41025e+01,  2.74999e+00,  2.47433e+00, ...,  8.17272e+01,
               -1.10226e+01, -3.96610e+01],
              ...,
              [ 2.31368e+01,  3.04407e+00,  3.61227e+00, ...,  2.67658e+01,
               -1.08499e+01, -7.09202e+00],
              [ 5.39127e+01,  3.17796e+00,  1.79006e+00, ..., -1.81374e+01,
               -1.65917e+01, -1.32651e+01],
              [ 2.67839e+01,  2.70001e+00,  3.35098e+00, ...,  1.36421e+01,
               -1.50580e+01, -7.83790e+00]]], shape=(4, 10000, 10))

    • x_c
      (chain, draw, x_c_dim_0)
      float64
      39.0 2.558 2.168 ... -9.94 -4.96

      array([[[  38.9981 ,    2.55818,    2.1681 , ...,   98.6768 ,
               1284.84   ,   -5.70191],
              [   9.56211,    2.66252,    2.3822 , ...,  -21.5883 ,
                 -5.27716,   -4.18692],
              [  19.7122 ,    2.36362,    3.10586, ...,  -96.2367 ,
                -35.4909 ,   -3.91073],
              ...,
              [  36.9361 ,    2.9664 ,    2.61065, ...,  -10.7237 ,
                 -4.34783,   -3.71794],
              [ -38.4519 ,    2.33649,    3.745  , ...,  -36.7467 ,
                 -4.9224 ,   -4.62733],
              [  19.829  ,    2.33427,    2.36371, ...,   88.8369 ,
                 -5.3189 ,   -2.79425]],
      
             [[  19.1122 ,    2.77597,    2.74671, ...,  -38.7813 ,
                -11.844  ,  -23.5007 ],
              [ 170.25   ,    2.6279 ,    2.39545, ..., -114.741  ,
                -10.7947 ,   -6.25547],
              [  22.1843 ,    2.68234,    2.31553, ...,  -24.5305 ,
                 -5.46523,   -2.95849],
      ...
              [ -73.1197 ,    2.68006,    2.20446, ...,  -17.6942 ,
                 -5.73777,   -4.20435],
              [  13.34   ,    2.50971,    2.49768, ...,  -39.5773 ,
                -10.1497 ,   -6.25722],
              [ -24.967  ,    2.48204,    2.06507, ..., -127.489  ,
                -13.5136 ,  -11.8541 ]],
      
             [[  13.8944 ,    2.51819,    2.33895, ..., 1119.25   ,
                -14.7083 ,   -4.42616],
              [   5.77262,    2.49245,    2.6532 , ...,  -26.9132 ,
                 -9.77818,   -3.57684],
              [  45.677  ,    2.39663,    2.23881, ...,   58.2514 ,
                 -7.07408,  -27.044  ],
              ...,
              [  16.7558 ,    2.59239,    3.11925, ...,   19.3374 ,
                 -7.01945,   -4.43243],
              [  38.5961 ,    2.71457,    1.76461, ...,  -12.2655 ,
                -10.9434 ,   -8.62686],
              [  19.3216 ,    2.36918,    2.94345, ...,   10.1957 ,
                 -9.93986,   -4.95961]]], shape=(4, 10000, 10))

    • y_c
      (chain, draw, y_c_dim_0)
      float64
      39.07 2.491 1.931 ... -10.33 -5.066

      array([[[  39.0666 ,    2.49112,    1.93134, ...,   99.4206 ,
               1302.76   ,   -5.85138],
              [   9.71507,    2.63405,    2.15352, ...,  -21.8346 ,
                 -5.44191,   -4.51193],
              [  19.5906 ,    2.24511,    2.85444, ...,  -97.0435 ,
                -37.0451 ,   -4.0918 ],
              ...,
              [  36.2472 ,    2.90341,    2.1873 , ...,  -10.726  ,
                 -4.44965,   -3.78589],
              [ -39.2476 ,    2.25144,    3.41016, ...,  -36.8423 ,
                 -5.0007 ,   -4.66264],
              [  20.5852 ,    2.29305,    2.15959, ...,   90.5002 ,
                 -5.34897,   -2.82634]],
      
             [[  19.5498 ,    2.71352,    2.54359, ...,  -38.6612 ,
                -12.2402 ,  -24.017  ],
              [ 176.577  ,    2.56361,    2.17534, ..., -113.413  ,
                -11.0628 ,   -6.08528],
              [  22.4601 ,    2.62666,    2.12662, ...,  -25.6929 ,
                 -5.76042,   -3.27536],
      ...
              [ -74.1182 ,    2.61571,    1.98335, ...,  -17.7185 ,
                 -5.97157,   -4.56864],
              [  13.4802 ,    2.42069,    2.30714, ...,  -39.892  ,
                -10.4061 ,   -6.18762],
              [ -25.2831 ,    2.44789,    1.80545, ..., -127.751  ,
                -14.092  ,  -12.2826 ]],
      
             [[  13.8206 ,    2.45003,    2.15033, ..., 1126.96   ,
                -15.5193 ,   -4.60126],
              [   5.91836,    2.41394,    2.3899 , ...,  -27.0929 ,
                 -9.79761,   -3.70761],
              [  45.9678 ,    2.34409,    2.07628, ...,   58.3331 ,
                 -7.49204,  -28.0244 ],
              ...,
              [  16.9356 ,    2.56105,    2.83576, ...,   19.5054 ,
                 -7.29924,   -4.53395],
              [  38.6274 ,    2.63856,    1.58322, ...,  -12.3748 ,
                -11.5006 ,   -9.09872],
              [  19.5264 ,    2.3027 ,    2.64772, ...,   10.1132 ,
                -10.3331 ,   -5.06614]]], shape=(4, 10000, 10))

  • Indexes: (9)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             9990, 9991, 9992, 9993, 9994, 9995, 9996, 9997, 9998, 9999],
            dtype='int64', name='draw', length=10000))

    • k_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='k_dim_0'))

    • dist_to_edge_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='dist_to_edge_dim_0'))

    • int_real_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='int_real_dim_0'))

    • phi_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='phi_dim_0'))

    • R_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='R_dim_0'))

    • x_c_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='x_c_dim_0'))

    • y_c_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='y_c_dim_0'))

  • Attributes: (4)

    created_at :

    2025-08-12T18:23:50.607928+00:00

    arviz_version :

    0.21.0

    inference_library :

    cmdstanpy

    inference_library_version :

    1.2.5

  
  
• sample_stats

  <xarray.Dataset> Size: 2MB
  Dimensions:          (chain: 4, draw: 10000)
  Coordinates:
    * chain            (chain) int64 32B 0 1 2 3
    * draw             (draw) int64 80kB 0 1 2 3 4 5 ... 9995 9996 9997 9998 9999
  Data variables:
      lp               (chain, draw) float64 320kB 455.8 451.5 ... 455.9 441.6
      acceptance_rate  (chain, draw) float64 320kB 0.9691 0.941 ... 0.9784 0.3328
      step_size        (chain, draw) float64 320kB 0.1824 0.1824 ... 0.1846 0.1846
      tree_depth       (chain, draw) int64 320kB 4 5 5 5 5 5 5 5 ... 4 5 5 4 5 5 5
      n_steps          (chain, draw) int64 320kB 31 31 31 31 31 ... 31 31 31 31 31
      diverging        (chain, draw) bool 40kB False False False ... False False
      energy           (chain, draw) float64 320kB -438.0 -442.5 ... -439.5 -431.0
  Attributes:
      created_at:                 2025-08-12T18:23:50.613913+00:00
      arviz_version:              0.21.0
      inference_library:          cmdstanpy
      inference_library_version:  1.2.5

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 10000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 9996 9997 9998 9999

      array([   0,    1,    2, ..., 9997, 9998, 9999], shape=(10000,))

  • Data variables: (7)
    • lp
      (chain, draw)
      float64
      455.8 451.5 450.3 ... 455.9 441.6

      array([[455.805, 451.478, 450.303, ..., 443.278, 445.157, 450.04 ],
             [448.474, 448.186, 448.89 , ..., 448.447, 448.657, 456.951],
             [448.382, 454.939, 448.659, ..., 454.089, 454.14 , 458.505],
             [449.564, 443.727, 446.883, ..., 452.326, 455.897, 441.568]],
            shape=(4, 10000))

    • acceptance_rate
      (chain, draw)
      float64
      0.9691 0.941 ... 0.9784 0.3328

      array([[0.969051, 0.941032, 0.824237, ..., 0.909056, 0.704401, 0.922582],
             [0.976227, 1.      , 0.829086, ..., 0.744463, 0.785643, 0.930751],
             [0.972075, 0.940094, 0.87217 , ..., 0.922485, 0.945976, 0.966526],
             [0.845398, 0.756874, 0.928666, ..., 0.855886, 0.97843 , 0.332781]],
            shape=(4, 10000))

    • step_size
      (chain, draw)
      float64
      0.1824 0.1824 ... 0.1846 0.1846

      array([[0.182397, 0.182397, 0.182397, ..., 0.182397, 0.182397, 0.182397],
             [0.173173, 0.173173, 0.173173, ..., 0.173173, 0.173173, 0.173173],
             [0.146456, 0.146456, 0.146456, ..., 0.146456, 0.146456, 0.146456],
             [0.184594, 0.184594, 0.184594, ..., 0.184594, 0.184594, 0.184594]],
            shape=(4, 10000))

    • tree_depth
      (chain, draw)
      int64
      4 5 5 5 5 5 5 5 ... 5 4 5 5 4 5 5 5

      array([[4, 5, 5, ..., 5, 4, 5],
             [5, 5, 5, ..., 5, 5, 5],
             [5, 5, 5, ..., 5, 5, 5],
             [4, 5, 5, ..., 5, 5, 5]], shape=(4, 10000))

    • n_steps
      (chain, draw)
      int64
      31 31 31 31 31 ... 31 31 31 31 31

      array([[31, 31, 31, ..., 31, 31, 31],
             [31, 31, 31, ..., 31, 31, 31],
             [31, 31, 31, ..., 31, 31, 31],
             [15, 31, 31, ..., 31, 31, 31]], shape=(4, 10000))

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 10000))

    • energy
      (chain, draw)
      float64
      -438.0 -442.5 ... -439.5 -431.0

      array([[-438.017, -442.458, -436.812, ..., -427.45 , -423.89 , -431.633],
             [-427.897, -432.709, -426.636, ..., -422.201, -436.194, -436.847],
             [-436.593, -431.605, -433.958, ..., -433.46 , -437.903, -438.896],
             [-431.851, -429.798, -427.946, ..., -433.911, -439.546, -430.961]],
            shape=(4, 10000))

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             9990, 9991, 9992, 9993, 9994, 9995, 9996, 9997, 9998, 9999],
            dtype='int64', name='draw', length=10000))

  • Attributes: (4)

    created_at :

    2025-08-12T18:23:50.613913+00:00

    arviz_version :

    0.21.0

    inference_library :

    cmdstanpy

    inference_library_version :

    1.2.5

  
  

The raw data from this object isn’t particularly helpful. Fortunately, we have a helper function in the extract module called extract_values(). You can extract the following parameters:

int_site : predicted paleointensity

sd_site : predicted standard deviation from best fit line of paleointensity

k : predicted curvature for each specimen

dist_to_edge: predicted \(D\) for each specimen (see Cych et al. 2021)

c: predicted line slope (relationship between intensity and curvature) for the site

phi : predicted \(\phi\) for each specimen (see Cych et al 2021.)

R : predicted \(R\) same as (\(1/\vec{k}\)) for each specimen

x_c : precited \(x_{c}\) for each specimen (see Cych et al 2021.)

y_c : predicted \(y_{c}\) for each specimen (see Cych et al 2021.)

Let’s extract int_site and plot a histogram. It should look the same as the result from site.histplot()

from BiCEP.extract import extract_values
pint_dist=extract_values(site.fit,'int_site')
plt.hist(pint_dist,bins=100,density=True);
plt.ylabel('Probability Density')
plt.xlabel('Paleointensity ($\mu$T)');

Getting Sampler Diagnostics

The \(\hat{R}\) sampler diagnostics can be obtained using the site.get_specimen_rhats() method. This gives the worst \(\hat{R}\) value from all of the specimens.

To check this for each specimen, we can check the specimen.rhat attribute.

print('Worst r_hat:',site.get_specimen_rhats())
print('hw126a2 r_hat:', specimen.rhat)
print('hw126a7 r_hat:', hw126a7.rhat)

Worst r_hat: 1.0004317552003188
hw126a2 r_hat: 1.0001539525687053
hw126a7 r_hat: 1.000453824652979

On this run, hw126a7 was the source of the largest \(\hat{R}\) - even though this is the worst, we find $ 1.0 < \hat{R} < 1.1 $ - which indicates that there were no problems with the sampling (may change between runs).

Choosing models and MCMC sampling options

Sometimes, high \(\hat{R}\) can be caused by sampling issues that can be fixed by using a different parameterization of the model. This is an optional parameter in site.BiCEP_fit(). The two models are BiCEP.model_circle_fast and BiCEP.model_circle_slow (found in the root BiCEP package). Default behaviour is to use the slow model when the number of specimens is less than 7 and the fast model when the number of specimens is greater than 7. The slow model can give more detailed results, but sometimes may perform worse. Try the other if the default doesn’t seem to be working. Other parameters for BiCEP_fit include n_samples which controls the number of samples (default 30000). Upping this number will take longer but improve the fit.

Using BiCEP to group by sample instead of site

This can be done by calling the data.switch_grouping() method. It doesn’t work on an individual sample/site basis (but this can be edited in the csv files BiCEP uses).