This page was generated from docs/Examples/Error_propagation/Cpx_Liq_Thermobarometry_Error_prop.ipynb. Interactive online version: Binder badge.

Python Notebook Download

Error propagation - Cpx-Liquid.

You need to install Thermobar once on your machine, if you haven’t done this yet, uncomment the line below (remove the #)

[2]:
#!pip install Thermobar
[3]:
# Import other python stuff
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import Thermobar as pt
pd.options.display.max_columns = None

Importing data

[4]:
out=pt.import_excel('Cpx_Liq_Error_prop_Feig2010_example.xlsx', sheet_name="Sheet1")
my_input=out['my_input']
myCpxs1=out['Cpxs']
myLiquids1=out['Liqs']

[5]:
## Check input looks reasonable
display(myCpxs1.head())
display(myLiquids1.head())
SiO2_Cpx TiO2_Cpx Al2O3_Cpx FeOt_Cpx MnO_Cpx MgO_Cpx CaO_Cpx Na2O_Cpx K2O_Cpx Cr2O3_Cpx Sample_ID_Cpx
0 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.319704 0 0.252186 0
1 51.578350 0.356860 3.600200 6.390180 0.171190 16.402000 20.388060 0.332720 0 0.186270 1
2 51.567875 0.256271 4.763807 4.229021 0.000000 16.969344 21.061352 0.316506 0 0.528929 2
3 52.083950 0.325983 3.919567 5.285317 0.168583 16.604683 20.364050 0.318067 0 0.377937 3
4 52.331846 0.406175 3.974610 6.039907 0.173175 16.297102 19.784783 0.388863 0 0.285975 4
SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq MnO_Liq MgO_Liq CaO_Liq Na2O_Liq K2O_Liq Cr2O3_Liq P2O5_Liq H2O_Liq Fe3Fet_Liq NiO_Liq CoO_Liq CO2_Liq Sample_ID_Liq
0 50.973845 0.489282 19.349173 5.327336 0 4.511373 9.161873 4.205536 0 0 0 5.12 0.0 0 0 0 0
1 53.642629 0.617371 19.320971 4.884986 0 3.257343 6.824314 5.057043 0 0 0 5.18 0.0 0 0 0 1
2 49.629083 0.365787 19.103075 5.302860 0 6.403377 11.615943 3.275783 0 0 0 5.05 0.0 0 0 0 2
3 51.624314 0.443871 18.451100 6.303071 0 6.132229 10.365057 3.685714 0 0 0 2.67 0.0 0 0 0 3
4 53.980586 0.817471 17.454986 6.744429 0 5.149100 9.030871 4.179543 0 0 0 2.25 0.0 0 0 0 4

Example 1: Uncertainty in a single input parameter for Cpx-only and Cpx-Liq barometry

  • Here, we consider the effect of adding just 5% noise to measured Na2O contents of Cpx (a fairly typical uncertainty resulting from EPMA analyses)

  • Because our liquids and cpx dataframes need to be the same size to feed into the calculate_Cpx_Liq functions, we also use the add_noise_sample_1phase to generate a dataframe of Liq compositions, however, we simply state noise_percent=0 so all the rows for each liquid are identical

[6]:
# Make liquid dataframe of the right size, no noise added
Liquids_only_noNoise=pt.add_noise_sample_1phase(phase_comp=myLiquids1,
                        noise_percent=0, duplicates=1000, err_dist="normal")

All negative numbers replaced with zeros. If you wish to keep these, set positive=False
[7]:
# Make Cpx dataframe, with 5% random noise added to Na2O
Cpx_5Na2O=pt.add_noise_sample_1phase(phase_comp=myCpxs1, variable="Na2O",
variable_err=5, variable_err_type="Perc", duplicates=1000, err_dist="normal")

Cpx_5Na2O.head()
All negative numbers replaced with zeros. If you wish to keep these, set positive=False
[7]:
SiO2_Cpx TiO2_Cpx Al2O3_Cpx FeOt_Cpx MnO_Cpx MgO_Cpx CaO_Cpx Na2O_Cpx K2O_Cpx Cr2O3_Cpx Sample_ID_Cpx_Num Sample_ID_Cpx
0 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.339825 0.0 0.252186 0.0 0.0
1 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.309925 0.0 0.252186 0.0 0.0
2 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.285936 0.0 0.252186 0.0 0.0
3 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.317906 0.0 0.252186 0.0 0.0
4 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.341860 0.0 0.252186 0.0 0.0
[8]:
# Plot what the Na content looks like
plt.hist(Cpx_5Na2O['Na2O_Cpx'].loc[Cpx_5Na2O['Sample_ID_Cpx_Num']==0], ec='k')
plt.xlabel('Na2O Cpx #0')
[8]:
Text(0.5, 0, 'Na2O Cpx #0')
../../_images/Examples_Error_propagation_Cpx_Liq_Thermobarometry_Error_prop_10_1.png

Lets first calculate Cpx-only pressures for these synthetic noisy Cpxs

[9]:
Out_5_noise_cpx_only=pt.calculate_cpx_only_press_temp(cpx_comps=Cpx_5Na2O,
                                    equationP="P_Put2008_eq32a",
                                        equationT="T_Put2008_eq32d", return_input=True)
Out_5_noise_cpx_only
[9]:
P_kbar_calc T_K_calc Delta_P_kbar_Iter Delta_T_K_Iter SiO2_Cpx TiO2_Cpx Al2O3_Cpx FeOt_Cpx MnO_Cpx MgO_Cpx CaO_Cpx Na2O_Cpx K2O_Cpx Cr2O3_Cpx Sample_ID_Cpx_Num Sample_ID_Cpx
0 7.423990 1504.919136 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.339825 0.0 0.252186 0.0 0.0
1 7.220311 1503.672762 4.547474e-13 3.865352e-12 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.309925 0.0 0.252186 0.0 0.0
2 7.056752 1502.665863 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.285936 0.0 0.252186 0.0 0.0
3 7.274696 1504.006386 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.317906 0.0 0.252186 0.0 0.0
4 7.437846 1505.003614 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.341860 0.0 0.252186 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
6995 7.288158 1510.951340 0.000000e+00 0.000000e+00 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.337092 0.0 0.377700 6.0 6.0
6996 7.292967 1510.984335 4.547474e-13 3.865352e-12 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.337779 0.0 0.377700 6.0 6.0
6997 7.165735 1510.110675 0.000000e+00 0.000000e+00 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.319615 0.0 0.377700 6.0 6.0
6998 7.341418 1511.316583 9.094947e-13 7.503331e-12 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.344694 0.0 0.377700 6.0 6.0
6999 7.146571 1509.978935 4.547474e-13 3.865352e-12 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.316879 0.0 0.377700 6.0 6.0

7000 rows × 16 columns

Lets plot the distribution of pressures for each Cpx

[10]:
fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(10, 10))

ax1.hist(Out_5_noise_cpx_only.loc[Out_5_noise_cpx_only['Sample_ID_Cpx_Num']==0, "P_kbar_calc"], bins=50, ec='k', density = True)
ax1.annotate("Cpx 1", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax2.hist(Out_5_noise_cpx_only.loc[Out_5_noise_cpx_only['Sample_ID_Cpx_Num']==1, "P_kbar_calc"], bins=50, ec='k', density = True)
ax2.annotate("Cpx 2", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax3.hist(Out_5_noise_cpx_only.loc[Out_5_noise_cpx_only['Sample_ID_Cpx_Num']==2, "P_kbar_calc"], bins=50, ec='k', density = True)
ax3.annotate("Cpx 3", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax4.hist(Out_5_noise_cpx_only.loc[Out_5_noise_cpx_only['Sample_ID_Cpx_Num']==3, "P_kbar_calc"], bins=50, ec='k', density = True)
ax4.annotate("Cpx 4", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax5.hist(Out_5_noise_cpx_only.loc[Out_5_noise_cpx_only['Sample_ID_Cpx_Num']==4, "P_kbar_calc"], bins=50, ec='k', density = True)
ax5.annotate("Cpx 5", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax6.hist(Out_5_noise_cpx_only.loc[Out_5_noise_cpx_only['Sample_ID_Cpx_Num']==5, "P_kbar_calc"], bins=50, ec='k', density = True)
ax6.annotate("Cpx 6", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)

ax1.set_xlabel('Pressure (kbar)')
ax2.set_xlabel('Pressure (kbar)')
ax4.set_xlabel('Pressure (kbar)')
ax3.set_xlabel('Pressure (kbar)')
ax6.set_xlabel('Pressure (kbar)')
ax5.set_xlabel('Pressure (kbar)')
ax1.set_ylabel('Probability Density')
ax3.set_ylabel('Probability Density')
ax5.set_ylabel('Probability Density')
plt.rcParams['figure.dpi']= 300
../../_images/Examples_Error_propagation_Cpx_Liq_Thermobarometry_Error_prop_14_0.png

We can also calculate statistics…

  • The first input is the thing you want to average (e.g. the calculated temperature in this instance), the second input is the thing you want to do the averaging based on (e.g. here average all rows with the same value of “Sample_ID_Liq_Num”

[11]:
Stats_T_K_cpx_only_5=pt.av_noise_samples_series(calc=Out_5_noise_cpx_only['T_K_calc'],
                                     sampleID=Out_5_noise_cpx_only['Sample_ID_Cpx_Num'])
Stats_T_K_cpx_only_5
[11]:
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 0.0 1000 1504.065074 1504.059012 0.665892 1506.150080 1501.700050
1 1.0 1000 1482.827485 1482.841381 0.596459 1484.663181 1480.767600
2 2.0 1000 1513.498688 1513.488535 0.604406 1515.522788 1511.592925
3 3.0 1000 1505.883677 1505.912476 0.694376 1508.855405 1503.531102
4 4.0 1000 1512.164707 1512.135692 0.926892 1515.248095 1509.528814
5 5.0 1000 1542.908695 1542.904504 1.236274 1546.890364 1539.202689
6 6.0 1000 1510.256541 1510.247912 0.762812 1512.385595 1507.692542
[47]:
Stats_P_kbar_cpx_only_5=pt.av_noise_samples_series(calc=Out_5_noise_cpx_only['P_kbar_calc'],
                                     sampleID=Out_5_noise_cpx_only['Sample_ID_Cpx_Num'])
Stats_P_kbar_cpx_only_5
[47]:
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 0.0 1000 7.288798 7.290636 0.108776 7.624785 6.864023
1 1.0 1000 5.692855 5.691012 0.106459 6.028094 5.330062
2 2.0 1000 7.377246 7.378669 0.107135 7.713221 6.993251
3 3.0 1000 7.089808 7.091796 0.108895 7.575986 6.775189
4 4.0 1000 8.383338 8.384756 0.140381 8.863947 7.912020
5 5.0 1000 10.872718 10.871257 0.184917 11.459944 10.224415
6 6.0 1000 7.182663 7.185273 0.112399 7.506819 6.872825

Lets do the same but for Cpx-Liq now. Remember, only Cpx has noise added

  • here, T=equation 33 from Putirka 2008, P=Equation 31 from Putirka 2008

[48]:
Out_5_noise_cpx_liq=pt.calculate_cpx_liq_press_temp(liq_comps=Liquids_only_noNoise, cpx_comps=Cpx_5Na2O,
                                    equationP="P_Put2008_eq31", equationT="T_Put2008_eq33", eq_tests=True)
Using Fe3FeT from input file to calculate Kd Fe-Mg

Each histogram shows the pressure distribution from a single Cpx-Liquid pair resulting from adding 5% error to just Na in Cpx.

  • Here we use the loc function to find all the rows where Cpx_Num==0, e.g. this is the first liquid in your excel spreadsheet.

  • If you wanted the 6th liquid instead, change ==0 to ==5 (as python starts its numbering from 0)

[49]:
fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(10, 10))

ax1.hist(Out_5_noise_cpx_liq.loc[Out_5_noise_cpx_liq['Sample_ID_Cpx_Num']==0, "P_kbar_calc"], bins=50, ec='k', density = True)
ax1.annotate("Liquid 1-Cpx 1", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax2.hist(Out_5_noise_cpx_liq.loc[Out_5_noise_cpx_liq['Sample_ID_Cpx_Num']==1, "P_kbar_calc"], bins=50, ec='k', density = True)
ax2.annotate("Liquid 2-Cpx 2", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax3.hist(Out_5_noise_cpx_liq.loc[Out_5_noise_cpx_liq['Sample_ID_Cpx_Num']==2, "P_kbar_calc"], bins=50, ec='k', density = True)
ax3.annotate("Liquid 3-Cpx 3", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax4.hist(Out_5_noise_cpx_liq.loc[Out_5_noise_cpx_liq['Sample_ID_Cpx_Num']==3, "P_kbar_calc"], bins=50, ec='k', density = True)
ax4.annotate("Liquid 4-Cpx 4", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax5.hist(Out_5_noise_cpx_liq.loc[Out_5_noise_cpx_liq['Sample_ID_Cpx_Num']==4, "P_kbar_calc"], bins=50, ec='k', density = True)
ax5.annotate("Liquid 5-Cpx 5", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)
ax6.hist(Out_5_noise_cpx_liq.loc[Out_5_noise_cpx_liq['Sample_ID_Cpx_Num']==5, "P_kbar_calc"], bins=50, ec='k', density = True)
ax6.annotate("Liquid 6-Cpx 6", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)

ax6.set_xlabel('Pressure (kbar)')
ax5.set_xlabel('Pressure (kbar)')
ax1.set_ylabel('Probability Density')
plt.rcParams['figure.dpi']= 300
../../_images/Examples_Error_propagation_Cpx_Liq_Thermobarometry_Error_prop_21_0.png

Calculate Statistics

[50]:
Stats_T_K_Cpx_Liq_5Na=pt.av_noise_samples_series(calc=Out_5_noise_cpx_liq['T_K_calc'],
                                     sampleID=Out_5_noise_cpx_liq['Sample_ID_Liq_Num'])
Stats_T_K_Cpx_Liq_5Na

[50]:
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 0.0 1000 1344.972079 1345.066097 2.375832 1351.859302 1334.858858
1 1.0 1000 1300.540818 1300.553805 2.177376 1306.980952 1292.523439
2 2.0 1000 1384.041287 1384.133098 2.586362 1391.667645 1374.029512
3 3.0 1000 1397.446766 1397.555216 2.659109 1408.299400 1389.294034
4 4.0 1000 1402.342547 1402.431341 2.676234 1410.894436 1392.658888
5 5.0 1000 1473.001786 1473.044707 2.988632 1481.914832 1461.717467
6 6.0 1000 1480.737532 1480.881133 3.166913 1489.359288 1471.540837
[51]:
Stats_P_kbar_Cpx_Liq_5Na=pt.av_noise_samples_series(calc=Out_5_noise_cpx_liq['P_kbar_calc'],
                                     sampleID=Out_5_noise_cpx_liq['Sample_ID_Liq_Num'])
Stats_P_kbar_Cpx_Liq_5Na
[51]:
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 0.0 1000 3.937346 3.946886 0.240857 4.635528 2.911816
1 1.0 1000 3.756044 3.757380 0.229345 4.434698 2.911861
2 2.0 1000 4.353471 4.362567 0.256145 5.108783 3.361773
3 3.0 1000 4.498447 4.509058 0.259736 5.558362 3.701910
4 4.0 1000 5.901296 5.909966 0.261498 6.737050 4.955068
5 5.0 1000 8.103374 8.107373 0.280832 8.941269 7.043299
6 6.0 1000 6.213997 6.227178 0.288307 6.998148 5.375994

Example 2 - published absolute 1 sigma values for all oxides

  • Here, we use the published 1 sigma values for each oxide in both the glass and cpx reported by Feig et al. 2010

  • We make a noisy dataframe of the same length for liquids and cpxs, then combine them into the function for iterating P and T

First, load in the published 1 sigma values

[52]:
out_err=pt.import_excel_errors('Cpx_Liq_Error_prop_Feig2010_example.xlsx', sheet_name="Sheet1")
myLiquids1_err=out_err['Liqs_Err']
myCpxs1_err=out_err['Cpxs_Err']
myinput_Out=out_err['my_input_Err']

## Check input looks reasonable
display(myCpxs1_err.head())
display(myLiquids1_err.head())

SiO2_Cpx_Err TiO2_Cpx_Err Al2O3_Cpx_Err FeOt_Cpx_Err MnO_Cpx_Err MgO_Cpx_Err CaO_Cpx_Err Na2O_Cpx_Err K2O_Cpx_Err Cr2O3_Cpx_Err P_kbar_Err T_K_Err Sample_ID_Cpx_Err
0 0.400044 0.032172 0.510460 0.357680 0.023249 0.538039 0.645164 0.049138 0 0.034646 0 0 0
1 0.628064 0.067227 0.622601 0.657476 0.050037 0.732816 1.050758 0.060948 0 0.062095 0 0 1
2 0.585596 0.032846 0.918768 0.506392 0.000000 0.562266 0.794310 0.107095 0 0.210114 0 0 2
3 0.165277 0.032709 0.240924 0.207266 0.051434 0.182666 0.355487 0.048979 0 0.058814 0 0 3
4 0.106509 0.040535 0.479962 0.428886 0.026673 0.431550 0.486264 0.056668 0 0.051372 0 0 4
SiO2_Liq_Err TiO2_Liq_Err Al2O3_Liq_Err FeOt_Liq_Err MnO_Liq_Err MgO_Liq_Err CaO_Liq_Err Na2O_Liq_Err K2O_Liq_Err Cr2O3_Liq_Err P2O5_Liq_Err H2O_Liq_Err Fe3Fet_Liq_Err NiO_Liq_Err CoO_Liq_Err CO2_Liq_Err Sample_ID_Liq_Err P_kbar_Err T_K_Err
0 0.326030 0.042448 0.221884 0.432684 0 0.207011 0.268045 0.213647 0 0 0 0.2560 0.0 0 0 0 0 0 0
1 0.328698 0.025182 0.242675 0.214755 0 0.163292 0.171717 0.085388 0 0 0 0.2590 0.0 0 0 0 1 0 0
2 0.476011 0.027083 0.235814 0.301696 0 0.174436 0.206030 0.229583 0 0 0 0.2525 0.0 0 0 0 2 0 0
3 0.442894 0.047517 0.234808 0.219780 0 0.171799 0.274484 0.259757 0 0 0 0.1335 0.0 0 0 0 3 0 0
4 0.445356 0.043963 0.316155 0.389311 0 0.134659 0.209045 0.453723 0 0 0 0.1125 0.0 0 0 0 4 0 0
[53]:
# Make noise based on reported absolute 1 sigma errors from Feig et al (2010) for liquids
Liquids_st_noise=pt.add_noise_sample_1phase(phase_comp=myLiquids1, phase_err=myLiquids1_err,
                                             phase_err_type="Abs", duplicates=1000, err_dist="normal")

# Make noise based on reported absolute 1 sigma errors from Feig et al (2010) for Cpx
Cpxs_st_noise=pt.add_noise_sample_1phase(phase_comp=myCpxs1, phase_err=myCpxs1_err,
                                             phase_err_type="Abs", duplicates=1000, err_dist="normal")

All negative numbers replaced with zeros. If you wish to keep these, set positive=False
All negative numbers replaced with zeros. If you wish to keep these, set positive=False
[54]:
Out_st_noise=pt.calculate_cpx_liq_press_temp(liq_comps=Liquids_st_noise, cpx_comps=Cpxs_st_noise,
                                    equationP="P_Put2008_eq31",
                                    equationT="T_Put2008_eq33", eq_tests=True)
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2115: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2115: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2115: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
Using Fe3FeT from input file to calculate Kd Fe-Mg
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2115: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
[55]:
Stats_T_K_CpxLiqPubErr=pt.av_noise_samples_series(Out_st_noise['T_K_calc'], Liquids_st_noise['Sample_ID_Liq_Num'])
Stats_T_K_CpxLiqPubErr
[55]:
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 0.0 1000 1344.282352 1344.310234 11.736914 1382.223983 1308.252768
1 1.0 1000 1298.645493 1299.964726 15.118184 1335.732684 1171.845161
2 2.0 1000 1381.243598 1383.698359 22.709192 1432.682914 1202.891430
3 3.0 1000 1396.811303 1397.093660 10.179831 1429.272558 1362.025647
4 4.0 1000 1401.665716 1401.300418 10.857033 1434.289423 1349.895560
5 5.0 1000 1473.147243 1472.916666 6.717644 1496.312115 1453.590474
6 6.0 1000 1480.494218 1480.649307 8.944667 1509.792188 1451.049517
[56]:
Stats_P_kbar_CpxLiqPubErr=pt.av_noise_samples_series(Out_st_noise['P_kbar_calc'], Liquids_st_noise['Sample_ID_Liq_Num'])
Stats_P_kbar_CpxLiqPubErr
[56]:
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 0.0 1000 3.878076 3.948044 0.870732 6.433104 0.779263
1 1.0 1000 3.577891 3.730529 1.354368 6.458861 -9.112942
2 2.0 1000 4.085652 4.373375 2.122981 8.397044 -14.144250
3 3.0 1000 4.453940 4.487837 0.866787 6.826525 0.854963
4 4.0 1000 5.810056 5.837714 0.906032 8.432087 1.650600
5 5.0 1000 8.102453 8.129936 0.474165 9.467956 6.756773
6 6.0 1000 6.232726 6.244151 0.509950 7.699720 4.696435

Recreate plot shown in manuscript for just 2 cpx-liq pairs from experiments

  • Here, on the LHS we plot cpx-liq pairs from the first row of the spreadsheet (ID=0)

  • On the RHS, we plot the pairs from row 4 (ID=3)

[57]:
fig, ((ax1, ax2)) = plt.subplots(1, 2, figsize=(12, 4))

ax1.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Liq_Num']==0, "P_kbar_calc"], bins=50, ec='k', density = True)
ax1.annotate("Liquid 1-Cpx 1", xy=(0.02, 0.95), xycoords="axes fraction", fontsize=9)
ax2.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Liq_Num']==3, "P_kbar_calc"], bins=50, ec='k', density = True)
ax2.annotate("Liquid 3-Cpx 3", xy=(0.02, 0.95), xycoords="axes fraction", fontsize=9)

#ax2.set_xlim([0, 6])
ax1.set_xlabel('Pressure (kbar)')
ax2.set_xlabel('Pressure (kbar)')
ax1.set_ylabel('Probability Density')
ax2.set_ylabel('Probability Density')
fig.savefig('Manuscript_2CpxNoises.png', dpi=300, transparent=True)
../../_images/Examples_Error_propagation_Cpx_Liq_Thermobarometry_Error_prop_33_0.png

Plot for 6 separate cpx-liq pairs

  • can average by _Cpx or _Liq, doesnt matter, as will match

[58]:
fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(10, 10))

ax1.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==0, "P_kbar_calc"], bins=50, ec='k', density = True)
ax1.annotate("Liquid 1-Cpx 1", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)


ax2.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==1, "P_kbar_calc"], bins=50, ec='k', density = True)
ax2.annotate("Liquid 2-Cpx 2", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)

ax3.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==2, "P_kbar_calc"], bins=50, ec='k', density = True)
ax3.annotate("Liquid 3-Cpx 3", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)

ax4.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==3, "P_kbar_calc"], bins=50, ec='k', density = True)
ax4.annotate("Liquid 4-Cpx 4", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)

ax5.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==4, "P_kbar_calc"], bins=50, ec='k', density = True)
ax5.annotate("Liquid 5-Cpx 5", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)

ax6.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==5, "P_kbar_calc"], bins=50, ec='k', density = True)
ax6.annotate("Liquid 6-Cpx 6", xy=(0.02, 0.9), xycoords="axes fraction", fontsize=9)

ax2.set_xlim([0, 6])
ax6.set_xlabel('Pressure (kbar)')
ax5.set_xlabel('Pressure (kbar)')
ax1.set_ylabel('Probability Density')
ax3.set_ylabel('Probability Density')
ax5.set_ylabel('Probability Density')
[58]:
Text(0, 0.5, 'Probability Density')
../../_images/Examples_Error_propagation_Cpx_Liq_Thermobarometry_Error_prop_35_1.png

Temperature distributions for first 6 pairs

[59]:
fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(10, 10))

ax1.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==0, "T_K_calc"], bins=50, ec='k', density = True)
ax2.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==1, "T_K_calc"], bins=50, ec='k', density = True)
ax3.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==2, "T_K_calc"], bins=50, ec='k', density = True)
ax4.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==3, "T_K_calc"], bins=50, ec='k', density = True)
ax5.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==4, "T_K_calc"], bins=50, ec='k', density = True)
ax6.hist(Out_st_noise.loc[Out_st_noise['Sample_ID_Cpx_Num']==5, "T_K_calc"], bins=50, ec='k', density = True)
ax6.set_xlabel('T (K)')
ax5.set_xlabel('T (K)')
ax1.set_ylabel('Probability Density')
ax3.set_ylabel('Probability Density')
ax5.set_ylabel('Probability Density')
[59]:
Text(0, 0.5, 'Probability Density')
../../_images/Examples_Error_propagation_Cpx_Liq_Thermobarometry_Error_prop_37_1.png

Plots with errorbars for each Cpx-Liq pair

[60]:
fig, (ax1) = plt.subplots(1, 1, figsize=(3, 2.5))

ax1.errorbar(Stats_T_K_CpxLiqPubErr['Mean_calc'],
             Stats_P_kbar_CpxLiqPubErr['Mean_calc'],
             xerr=Stats_T_K_CpxLiqPubErr['St_dev_calc'],
             yerr=Stats_P_kbar_CpxLiqPubErr['St_dev_calc'],
             fmt='d', ecolor='k', elinewidth=0.8, mfc='cyan', ms=5, mec='k')
ax1.set_ylabel('Calculated P (kbar)')
ax1.set_xlabel('Calculated T (K)')


[60]:
Text(0.5, 0, 'Calculated T (K)')
../../_images/Examples_Error_propagation_Cpx_Liq_Thermobarometry_Error_prop_39_1.png
[25]:
Stats_P_kbar_CpxLiqPubErr
[25]:
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 0.0 1000 3.875331 3.897954 0.838025 6.452778 0.464844
1 1.0 1000 3.570149 3.813306 1.640338 5.931074 -20.725732
2 2.0 1000 4.061145 4.403422 2.120289 8.365879 -11.312059
3 3.0 1000 4.458786 4.515581 0.840769 6.767782 1.166771
4 4.0 1000 5.886142 5.920240 0.912282 8.860778 1.436836
5 5.0 1000 8.091746 8.073837 0.460574 9.480674 6.397698
6 6.0 1000 6.226063 6.260299 0.501330 7.886687 4.558481
[ ]: