Help with plots

[Swan-chip]

Hello all

Real beginner with pyleoclim but really enjoying. Struggling with plots though.

1. Using the pyleoclim.core.psds.PSD.plot() is there a way of using linear on the x-axis bu log-scale on the Y-axis?

2. On a related question, is there a way of extracting the numeric values from the PSD.signif_test() method?

I can find the PSD.signif_qs attribute, which is a MultiplePSD object that displays the confidence level graphically when called, but I have found no way to get to the underlying values. I have found some of the formatting features associated with the plots frustratingly unpredictable (problem with my install?) so it would be great to be able to extract all the data to plot using more familiar modules that I can tailor more precisely.

Many thanks in advance :-)

[Swan-chip]

I have answered my own question regards (2). The values can be extracted from the axes object using

axes.get_children()[index].get_xdata()

Regards (1), I have not found a way of using mixed axes directly using pyleoclim, and would still be interested, but now I can get all the values out I can work around this.

[CommonClimate]

Hello @Swan-chip and welcome to Pyleoclim! These are good questions and I need to make a note of adding that to a tutorial at some point. It's always easier with some code to demonstrate, so let's start from the docstring from the PSD.signif_test() method:

import pyleoclim as pyleo
soi = pyleo.utils.load_dataset('SOI')
psd = soi.standardize().spectral('mtm',settings={'NW':2})
psd_sim = psd.signif_test(number=20)

1. Using the pyleoclim.core.psds.PSD.plot() is there a way of using linear on the x-axis bu log-scale on the Y-axis?

Assuming you want to use the frequency axis for this instead of period:

fig, ax = psd_sim.plot(in_loglog=False, in_period=False)
ax.set_yscale('log')

You can also hack the source code for the plot method for inspiration (though it is admittedly not for the faint of heart!).

2. On a related question, is there a way of extracting the numeric values from the PSD.signif_test() method?

You will find them in psd_sim.signif_qs.psd_list[0].amplitude

The key here is that the MultiplePSD object is a list of PSDs. One accesses its elements via psd_list, which has as many elements as the qs argument in signif_test() (only 1 by default, the 95% quantile).

Please let us know if you have other questions! Your solution of pulling the data from the figure axes is resourceful, but a tad more cumbersome. It's rare to have to dig down to those values, but if that becomes more commonplace we can find ways to better expose them. Note also that the PSD object does a few things besides plotting, and in that case you did not need to plot to access the values. If you were looping over many/long series, these differences would become meaningful.

[Swan-chip]

Thank you so much for this, really helpful.

I have a closely related question. I am trying to understand a failure to reproduce confidence levels on a time series between my code and that written on Matlab using the RedConf script. I would like to check each step in the process. Is there a way of outputting and plotting the fitted AR1 noise itself rather than a confidence level based on the chi-squared distribution of that noise?

Thank you!

[Swan-chip]

Thank you so much for this, really helpful.

I have a closely related question. I am trying to understand a failure to reproduce confidence levels on a time series between my code and that written on Matlab using the RedConf script. I would like to check each step in the process. Is there a way of outputting and plotting the fitted AR1 noise itself rather than a confidence level based on the chi-squared distribution of that noise?

Thank you!

[CommonClimate]

Hello @Swan-chip
I am not familiar with RedConf. Is it this one?

Pyleoclim has several ways of assessing the significance of spectral peaks. Since you are interested in AR(1) benchmarks, know that there are two ways:

• asymptotic limit (ar1asym): uses the theoretical AR(1) spectrum with chi^2 intervals
• direct simulation (ar1sim), which obtains 95% quantiles at each frequency from a set of simulated AR(1) spectra.

Details of the testing implementation are here. More background on using surrogates (both parametric and non-parametric) is here.

Question for you: where do you see the significance of spectral peaks being particularly consequential? We are currently testing another method to establish significance, and would like some real-world examples where the result may matter.

[Swan-chip]

Thank you for another helpful reply 😊

Yes it is the RedConf you found. I guess to get the underlying noise I can use:

pyleo.utils.tsmodel.ar1_sim

Which I could then presumably add into the PSD.signif_qs by some means.

You might wonder why I do not simply use qs = 0.5 and compare the pyleoclim with matlab confidence levels. However, I already know these diverge and I wonder whether the difference comes in the original fitting of noise or the calculation of chi-squared confidence levels. The pyleo asymptotic confidence levels are MUCH lower than the direct simulation (AR1sim). They place just about every peak >> 99% significance, which I think I can probably just discard as unhelpful.

Interestingly, the overall amplitudes of the two plots are quite different, despite having identical peaks and trough wavelengths and relative amplitudes. I wonder whether this is a clue, but I am clutching at straws. I have attached the pyleo and matlab plots of the same data purportedly with comparable analysis.

To answer your question, I don’t think I can give much useful feedback on the significance as I am only dabbling in this to test reproducibility. I am reviewing a manuscript for which there are other important aspects on which I do actually have some expertise! It has been fascinating, slight struggle but possibly relevant to postdoc I'm starting in May. Pyleoclim is excellently supported with tutorials and docs and of course this forum!

Any further thoughts on the comparison welcome. I am still stuck with a problem of WHY I am reproducing such lower significance with pyleo.

Many thanks

[Swan-chip]

So, have learnt a lot since writing the above, answering much of my previous message. Followed the code through back to statsmodels.tsa.arima_process.arma_generate_sample and now better understand what is going on. Really fun! Importantly, by recreating the process in pyleoclim, I now think I understand the discrpeancy between the pyloeclim and the Matlab analysis. Would be really interested to know whether anyone thinks the following is problematic? As I explain below, I am not sold on it.

The procedure is as follows:

1. Create multiple noise simulations using method = "ar1sim". Simulations are fit to time series using SurrogateSeries.from_series() method.

2. Create a mean of those simulations.

3. Calculate spectra, method = "mtm", settings={"NW":2}

4. Create a best fit polynomial from that mean.

5. Calculate 95% confidence levels from polynomial using a chi-squared distribution with 6 degrees of freedom:

# value = an amplitude value from mean noise spectra
# dof = 6 degrees of freedom

chi_crit_95 = scipy.stats.chi2.ppf(0.95, dof)

value_at_95_conf_level = value * (chi_crit_95/dof)

This reproduces matlab-like results and appears to follow the same workflow. However I have two concerns.

First, why calculate chi-squared estimates when we can just take the quantiles of the simulations, which incidentally do not agree with chi-squared calculated values in the above example?

Second, the result is sensitive to degrees of freedom. Traditionally, calculating chi-squared confidence levels from AR1 red noise associated with the MTM approach, the chi-squared distribution has 2 x tapers degrees of freedom, which I believe is 6 for NW=2. However, I don't think the simulation is based strictly on the MTM process (?) so is this number of degrees of freedom still valid?

In summary, would simply taking the percentil of the actual simulations be a simpler and more robust approach, as used in my previous post? Should the alternative described above be considered unsound.

Many thanks in advance for any thoughts and discussion!

[CommonClimate]

Hi! I am glad your process of tinkering is leading you to better understand what is happening. That's the way it should be.

It appears that your current workflow needlessly mixes the asymptotic approach and the simulation-based (Monte Carlo) approach. If you are not limited by computation, I recommend using ar1sim and pushing the number of simulations to >=1000, and simply getting the quantiles from that. Assuming your spectrum was stored in my_psd, this would look like:

psd_sig = my_psd.signif_test(method='ar1sim', qs = [.95, .99],number=1000)
psd_sig.plot()

[CommonClimate]

Indeed, your computation is sensitive to DOFs, which should depend on the autocorrelation -- your pseudocode quotes a fixed DOF = 6, which doesn't sound right.

[Swan-chip]

I am getting quite into this now. I cannot understand why the PSD.signif_test(method = "ar1asym") is not producing results that seem right or that I can reproduce elsewhere. I have been through the maths and do have a question or two but I think to keep this simple I have sent some code, which uses the built-in demonstration data. This doesn't look as bad as it does on my data but it doesn't look right to me. Very interested in your thoughts.

import pyleoclim as pyleo
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

lr04 = pyleo.utils.load_dataset('LR04')
lr04_s = lr04.slice([0,3000])

fig, ax = lr04_s.standardize().interp().spectral(method = "mtm", settings = {"NW":2}).signif_test(method = "ar1asym", qs = [0.95, 0.99]).plot()

Many thanks!

[CommonClimate]

I am not happy with the way the asymptotic spectra have been implemented ; it is ripe for an overhaul. Here is a comparison to the direct simulation (ar1sim):

the code is:

import pyleoclim as pyleo
import matplotlib.pyplot as plt
lr04 = pyleo.utils.load_dataset('LR04')
lr04_s = lr04.slice([0,3000])
psd = lr04_s.standardize().interp().spectral(method = "mtm", settings = {"NW":2})
qs = [0.95, 0.99]
psd_sig = list()
psd_sig.append(psd.signif_test(method = "ar1asym",qs = qs))
psd_sig.append(psd.signif_test(method = "ar1sim",number=5000, qs = qs))
methods = ['Asymptotic', "Monte Carlo"]

fig, axs = plt.subplots(2,1,sharex=True,figsize=(8,8))
axs=axs.flatten()
for i, psd in enumerate(psd_sig):
    psd.plot(ax=axs[i], title=methods[i])

I would recommend sticking to direct simulation if you can spare the CPU time ; will post updates when we have a better alternative. Note that we are also in the process of implementing an alternative that does not require interpolation, and has a number of other desirable statistical properties. I will focus on this in the next few weeks.

[CommonClimate]

We had a development meeting today and I confirm that ar1asym is not to be used for this. It only makes sense with cwt as the spectral method, which is only sensible on dashboard(). See #645