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 #)

[1]:
#!pip install Thermobar
[2]:
# 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
pt.__version__
[2]:
'1.0.30'

Importing data

[3]:
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']

[4]:
## 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 Feig2010_1
1 51.578350 0.356860 3.600200 6.390180 0.171190 16.402000 20.388060 0.332720 0 0.186270 Feig2010_2
2 51.567875 0.256271 4.763807 4.229021 0.000000 16.969344 21.061352 0.316506 0 0.528929 Feig2010_3
3 52.083950 0.325983 3.919567 5.285317 0.168583 16.604683 20.364050 0.318067 0 0.377937 Feig2010_4
4 52.331846 0.406175 3.974610 6.039907 0.173175 16.297102 19.784783 0.388863 0 0.285975 Feig2010_5
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 Feig2010_1
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 Feig2010_2
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 Feig2010_3
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 Feig2010_4
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 Feig2010_5

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

[5]:
# 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
[6]:
# 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
[6]:
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.339916 0 0.252186 0 Feig2010_1
1 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.336336 0 0.252186 0 Feig2010_1
2 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.310841 0 0.252186 0 Feig2010_1
3 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.350643 0 0.252186 0 Feig2010_1
4 51.075956 0.325086 4.892502 5.714957 0.1657 16.897042 20.466034 0.355440 0 0.252186 0 Feig2010_1
[7]:
# 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')
[7]:
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

[8]:
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
[8]:
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.424612 1504.922927 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.339916 0 0.252186 0 Feig2010_1
1 7.400235 1504.774210 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.336336 0 0.252186 0 Feig2010_1
2 7.226553 1503.711082 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.310841 0 0.252186 0 Feig2010_1
3 7.497637 1505.367706 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.350643 0 0.252186 0 Feig2010_1
4 7.530282 1505.566177 0.000000e+00 0.000000e+00 51.075956 0.325086 4.892502 5.714957 0.165700 16.897042 20.466034 0.355440 0 0.252186 0 Feig2010_1
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
6995 7.155413 1510.039721 4.547474e-13 3.865352e-12 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.318141 0 0.377700 6 Feig2010_7
6996 7.198844 1510.338189 4.547474e-13 3.865352e-12 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.324343 0 0.377700 6 Feig2010_7
6997 7.239387 1510.616622 4.547474e-13 3.865352e-12 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.330130 0 0.377700 6 Feig2010_7
6998 7.221420 1510.493252 4.547474e-13 3.865352e-12 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.327566 0 0.377700 6 Feig2010_7
6999 7.294841 1510.997189 4.547474e-13 3.637979e-12 52.498000 0.299767 3.517867 5.903067 0.178867 17.869733 19.342233 0.338046 0 0.377700 6 Feig2010_7

7000 rows × 16 columns

Lets plot the distribution of pressures for each Cpx

[9]:
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”

[10]:
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
[10]:
Sample # averaged Mean_calc Median_calc St_dev_calc St_dev_calc_IQR Max_calc Min_calc
0 0 1000 1504.088131 1504.063057 0.626372 0.640921 1505.989933 1501.773707
1 1 1000 1482.793286 1482.794223 0.616343 0.604748 1484.907167 1480.660467
2 2 1000 1513.541879 1513.536840 0.606123 0.620552 1515.532491 1511.818183
3 3 1000 1505.895229 1505.906350 0.700749 0.701374 1507.925028 1503.692797
4 4 1000 1512.124989 1512.126904 0.943300 0.917292 1514.870823 1508.698241
5 5 1000 1542.980662 1542.959607 1.260798 1.231232 1547.818880 1539.170773
6 6 1000 1510.267668 1510.241719 0.750436 0.743942 1512.552378 1507.687625
[ ]:
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

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

[ ]:
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)

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)

[ ]:
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

Calculate Statistics

[ ]:
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

[ ]:
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

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

[ ]:
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())

[ ]:
# 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")

[ ]:
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)
Out_st_noise['Sample_ID_Liq']=Liquids_st_noise['Sample_ID_Liq']
[ ]:
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
[ ]:
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

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)

[ ]:
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)

Plot for 6 separate cpx-liq pairs

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

[ ]:
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')

Temperature distributions for first 6 pairs

[ ]:
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')

Plots with conventional errorbars for each Cpx-Liq pair

  • Note, although this is typically how thermobarometric data is displayed, it isnt very realistic, our error propagation shows that error ellipses would be better - see the example below

[ ]:
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)')


Error ellipses instead

[ ]:
Out_st_noise.columns[0:100]
[ ]:
Out_st_noise['Sample_ID_Liq'].unique()[0]
[ ]:
fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 4))

# you have to plot something on the axis before you try to add the confidence intervals

ax0.plot(Stats_T_K_CpxLiqPubErr['Mean_calc'],
             Stats_P_kbar_CpxLiqPubErr['Mean_calc'], 'ok', mfc='red', ms=4)


for sam in Out_st_noise['Sample_ID_Liq'].unique():
    Calc_PT_16=Out_st_noise.loc[Out_st_noise['Sample_ID_Liq']==sam].reset_index(drop=True)

    y=Calc_PT_16['P_kbar_calc'].loc[Calc_PT_16['P_kbar_calc']>-10].values
    x=Calc_PT_16['T_K_calc'].loc[Calc_PT_16['P_kbar_calc']>-10].values
    pt.matplotlib_confidence_ellipse(x, y, ax0, n_std=1,
                    edgecolor='r', lw=0.5, alpha=0.8)



ax0.set_xlabel('T (K)')
ax0.set_ylabel('P (kbar)')

ax1.plot(Stats_T_K_CpxLiqPubErr['Mean_calc'],
             Stats_P_kbar_CpxLiqPubErr['Mean_calc'], 'ok', mfc='red', ms=4)


for sam in Out_st_noise['Sample_ID_Liq'].unique():

    Calc_PT_16=Out_st_noise.loc[Out_st_noise['Sample_ID_Liq']==sam].reset_index(drop=True)

    y=Calc_PT_16['P_kbar_calc'].loc[Calc_PT_16['P_kbar_calc']>-10].values
    x=Calc_PT_16['T_K_calc'].loc[Calc_PT_16['P_kbar_calc']>-10].values
    if sam==Out_st_noise['Sample_ID_Liq'].unique()[0]:
        pt.matplotlib_confidence_ellipse(x, y, ax1, n_std=1,
                    edgecolor='red', lw=0.5, alpha=0.8, label='1σ MC simulations')
    else:
        pt.matplotlib_confidence_ellipse(x, y, ax1, n_std=1,
                    edgecolor='red', lw=0.5, alpha=0.8)






ax1.set_xlabel('T (K)')
ax1.set_ylabel('P (kbar)')
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='o', ecolor='k', elinewidth=0.8, mfc='red', ms=5, mec='k', label='Conventional error bars')
ax1.legend()
ax0.set_title('Cpx-Liq Error ellipses')
ax1.set_title('Cpx-Liq Error bars vs. ellipses')

fig.savefig('confidence_intervals_vs_errorbars.png', dpi=200)
[ ]: