Investigation of double-peak separation

Contents

Investigation of double-peak separation#

%load_ext autoreload
%autoreload 2

import numpy as np
import arviz as az
import pymc as pm
from pathlib import Path
from peak_performance import pipeline as pl, models, plots
from matplotlib import pyplot as plt

Load and inspect raw intensity data#

timeseries = np.load(Path(r"..\example\A2t2R1Part1_132_85.9_86.1.npy"))

fig, ax = plt.subplots()
ax.scatter(timeseries[0], timeseries[1], marker="x", color="black")
ax.set(
    xlabel="time / h",
    ylabel="intensity / a.u.",
)
plt.show()

../_images/61067a3e04f13822acc32c624c02399685631d7efb9567108be7cad36fc8ec16.png

Define a peak model#

pmodel = models.define_model_double_normal(
    time=timeseries[0],
    intensity=timeseries[1]
)
pmodel.to_graphviz()

c:\Users\osthege\AppData\Local\mambaforge\envs\pepe\Lib\site-packages\pymc\data.py:273: FutureWarning: ConstantData is deprecated. All Data variables are now mutable. Use Data instead.
  warnings.warn(

../_images/376fc13af4a8a457acf54888535cb0041772d03da347693efb51db75aabae6d8.svg

Inspect priors for separation and retention time#

fig, axs = plt.subplots(dpi=200, ncols=2, figsize=(12, 6))
ax = axs[0]
ax.hist(pm.draw(pmodel["separation"], draws=50000), bins=100, density=True)
ax.set(ylabel="density / -", xlabel="$\mathrm{\Delta t\ /\ min}$", title="separation")

ax = axs[1]
ax.hist(pm.draw(pmodel["offset"], draws=50000), bins=100, density=True)
ax.set(ylabel="density / -", xlabel="$\mathrm{\Delta t\ /\ min}$", title="offset")
plt.show()

../_images/edf66c24522d4ac9a6fbe08c449080cdc7fdc7576583bb0383d3888ad69a7a02.png

Run MCMC sampling and inspect results#

idata = pl.sampling(pmodel, tune=6_000, draws=2000)
idata = pl.posterior_predictive_sampling(pmodel, idata)
idata

Sampling: [L, baseline_intercept, baseline_slope, height, meanmean, noise, separation, std]

100.00% [32000/32000 00:03<00:00 Sampling chains: 0, Divergences: 0]

Sampling: [L]

arviz.InferenceData

posterior

posterior_predictive
sample_stats

prior

prior_predictive
observed_data
constant_data

warmup_posterior

warmup_sample_stats

az.plot_trace(idata, var_names=["meanmean", "separation", "offset", "mean", "std"], combined=True)
plt.tight_layout()
plt.show()

../_images/1d44f7dacafb388af3d69486340cbdbd8de209834ca04882911586d4d6d4d291.png

summary = az.summary(idata, var_names=["~y", "~baseline", "offset"])
summary

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
area[0]	316.879	28.606	260.096	368.756	0.267	0.189	11447.0	5609.0	1.0
area[1]	674.658	26.019	628.417	725.552	0.247	0.175	11094.0	6189.0	1.0
baseline_intercept	1115.569	39.801	1041.960	1192.575	0.752	0.533	2805.0	3543.0	1.0
baseline_slope	-21.651	3.153	-27.610	-15.693	0.060	0.042	2801.0	3638.0	1.0
height[0]	774.967	64.153	658.831	899.934	0.748	0.529	7341.0	6375.0	1.0
height[1]	1762.363	64.785	1639.151	1884.930	0.731	0.517	7857.0	6440.0	1.0
height_log__[0]	6.649	0.083	6.491	6.803	0.001	0.001	7341.0	6375.0	1.0
height_log__[1]	7.474	0.037	7.407	7.546	0.000	0.000	7857.0	6440.0	1.0
mean[0]	11.731	0.015	11.703	11.760	0.000	0.000	5173.0	5319.0	1.0
mean[1]	12.433	0.007	12.422	12.446	0.000	0.000	12180.0	6652.0	1.0
meanmean	12.082	0.008	12.067	12.098	0.000	0.000	5707.0	5689.0	1.0
noise	118.698	8.156	104.081	134.331	0.086	0.061	9008.0	4975.0	1.0
noise_log__	4.774	0.068	4.649	4.903	0.001	0.001	9008.0	4975.0	1.0
offset[0]	-0.351	0.008	-0.366	-0.336	0.000	0.000	5682.0	5867.0	1.0
offset[1]	0.351	0.008	0.336	0.366	0.000	0.000	5682.0	5867.0	1.0
separation	0.702	0.016	0.672	0.732	0.000	0.000	5682.0	5867.0	1.0
separation_log__	-0.354	0.023	-0.396	-0.311	0.000	0.000	5682.0	5867.0	1.0
sn[0]	6.561	0.713	5.195	7.862	0.008	0.006	7691.0	5958.0	1.0
sn[1]	14.918	1.163	12.746	17.075	0.013	0.009	8553.0	5860.0	1.0
std[0]	0.164	0.017	0.134	0.199	0.000	0.000	6537.0	5798.0	1.0
std[1]	0.153	0.007	0.140	0.165	0.000	0.000	7141.0	6114.0	1.0
std_log__[0]	-1.814	0.106	-2.006	-1.612	0.001	0.001	6537.0	5798.0	1.0
std_log__[1]	-1.879	0.044	-1.963	-1.800	0.001	0.000	7141.0	6114.0	1.0

plots.plot_posterior_predictive(
    identifier="peak_fit",
    time=idata.constant_data.time.values,
    intensity=idata.constant_data.intensity.values,
    path=None,
    idata=idata,
    discarded=False,
)

../_images/1336e4c869d272220cf493dc5f89136a8f7e9f7d8cbae958318599b777317e53.png

az.plot_pair(idata, var_names="noise,height,std,mean,meanmean".split(","), kind="hexbin");

../_images/d62b33d4b11270808e263c3857bdda985315f5c50490406d7d4b9bbd317cb0cf.png

ax = az.plot_ppc(idata, kind="cumulative")
fig = plt.gcf()
ax.set(
    ylabel="cumulative density / -",
    xlabel="intensity / a.u.",
)
plt.tight_layout()

../_images/93524136d5752c6bd8434b80cac2695a0952491ceb1c17d816cb306eed3bfdd6.png

%load_ext watermark
%watermark -idu

Last updated: 2024-05-11T14:18:49.035894+02:00