Skip to content
Snippets Groups Projects
Commit 52325325 authored by Manuel Marschall's avatar Manuel Marschall
Browse files

update to paper quality plots

parent 27472ea6
Branches
No related tags found
No related merge requests found
...@@ -103,8 +103,8 @@ def main(): ...@@ -103,8 +103,8 @@ def main():
# contradiction case # contradiction case
# x = [0.78, 0.84] # measurement data --> acc rate: 2.13 # x = [0.78, 0.84] # measurement data --> acc rate: 2.13
for measurement, label in zip([[0.95, 1.05], [1.198, 1.203], [0.78, 0.84]], for measurement, label in zip([[0.95, 1.05], [1.198, 1.203], [0.78, 0.84]],
["ideal", "highly_informative", "contradiction"]): ["opt", "info", "contra"]):
experiment_run = experiment + f"case_{label}/" experiment_run = experiment
if not os.path.exists(experiment_run): if not os.path.exists(experiment_run):
os.makedirs(experiment_run) os.makedirs(experiment_run)
...@@ -116,8 +116,9 @@ def main(): ...@@ -116,8 +116,9 @@ def main():
print(f"actual acceptance rate: {len(rejection_result['Y_accept'])/n_samples*100:.4g}%") print(f"actual acceptance rate: {len(rejection_result['Y_accept'])/n_samples*100:.4g}%")
fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"], xlim=(0.7, 1.5), ylim=(0, 10)) fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"], xlim=(0.7, 1.5), ylim=(0, 10),
fig.savefig(experiment_run + "x_samples.png") fontsize=22)
fig.savefig(experiment_run + "x_" + f"case_{label}.pdf")
result_dict = { result_dict = {
"PRIOR": model["variables"]["Y"]["samples"], "PRIOR": model["variables"]["Y"]["samples"],
"REJECTION": rejection_result["Y_accept"], "REJECTION": rejection_result["Y_accept"],
...@@ -128,13 +129,13 @@ def main(): ...@@ -128,13 +129,13 @@ def main():
sampler = mcmc_reference(model, nwalker=32, nsamples=10000) sampler = mcmc_reference(model, nwalker=32, nsamples=10000)
fig = plot_mcmc_chains(sampler.get_chain(), labels=[key for key, value in model["variables"].items()]) fig = plot_mcmc_chains(sampler.get_chain(), labels=[key for key, value in model["variables"].items()])
fig.savefig(experiment_run + "mcmc.png") fig.savefig(experiment_run + "mcmc.pdf")
# Adjust the MCMC sampler chains to discard ~ 2xtau and thin by tau/2 # Adjust the MCMC sampler chains to discard ~ 2xtau and thin by tau/2
# tau = sampler.get_autocorr_time() # tau = sampler.get_autocorr_time()
# print(tau) # print(tau)
result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0] result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0]
fig = plot_result_pdfs(result_dict, xlim=(-0.7, 0.7)) fig = plot_result_pdfs(result_dict, xlim=(-0.7, 0.7))
fig.savefig(experiment_run + "result.png") fig.savefig(experiment_run + "result_" + f"case_{label}.pdf")
if __name__ == "__main__": if __name__ == "__main__":
......
...@@ -167,7 +167,7 @@ def generate(experiment: str, n_samples: int, bootstrap: int): ...@@ -167,7 +167,7 @@ def generate(experiment: str, n_samples: int, bootstrap: int):
res_dict["acc_rate"].append(len(rejection_result["Y_accept"])/n_samples*100) # type: ignore res_dict["acc_rate"].append(len(rejection_result["Y_accept"])/n_samples*100) # type: ignore
fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"]) fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"])
fig.savefig(experiment_run + "x_samples.png") fig.savefig(experiment_run + "x_samples.pdf")
result_dict = { result_dict = {
"PRIOR": model["variables"]["Y"]["samples"], "PRIOR": model["variables"]["Y"]["samples"],
"REJECTION": rejection_result["Y_accept"], "REJECTION": rejection_result["Y_accept"],
...@@ -186,7 +186,7 @@ def generate(experiment: str, n_samples: int, bootstrap: int): ...@@ -186,7 +186,7 @@ def generate(experiment: str, n_samples: int, bootstrap: int):
# result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0] # result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0]
fig = plot_result_pdfs(result_dict) fig = plot_result_pdfs(result_dict)
fig.savefig(experiment_run + "result.png") fig.savefig(experiment_run + "result.pdf")
print(f"posterior mean: {np.mean(rejection_result['Y_accept']):.4e}") print(f"posterior mean: {np.mean(rejection_result['Y_accept']):.4e}")
print(f"posterior uncertainty: {np.std(rejection_result['Y_accept'], ddof=1):.4e}") print(f"posterior uncertainty: {np.std(rejection_result['Y_accept'], ddof=1):.4e}")
......
...@@ -140,7 +140,7 @@ def main(exp: int = 1): ...@@ -140,7 +140,7 @@ def main(exp: int = 1):
print(f"actual acceptance rate: {len(rejection_result['Y_accept'])/n_samples*100:.4g}%") print(f"actual acceptance rate: {len(rejection_result['Y_accept'])/n_samples*100:.4g}%")
fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"], xlim=(0.9, 1.3) if exp == 2 else None) fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"], xlim=(0.9, 1.3) if exp == 2 else None)
fig.savefig(experiment + "x_samples.png") fig.savefig(experiment + "x_samples.pdf")
result_dict = { result_dict = {
"PRIOR": model["variables"]["Y"]["samples"], "PRIOR": model["variables"]["Y"]["samples"],
"REJECTION": rejection_result["Y_accept"], "REJECTION": rejection_result["Y_accept"],
...@@ -151,13 +151,13 @@ def main(exp: int = 1): ...@@ -151,13 +151,13 @@ def main(exp: int = 1):
sampler = mcmc_reference(model, nwalker=32, nsamples=50000) sampler = mcmc_reference(model, nwalker=32, nsamples=50000)
fig = plot_mcmc_chains(sampler.get_chain(), labels=[key for key, value in model["variables"].items()]) fig = plot_mcmc_chains(sampler.get_chain(), labels=[key for key, value in model["variables"].items()])
fig.savefig(experiment + "mcmc.png") fig.savefig(experiment + "mcmc.pdf")
# Adjust the MCMC sampler chains to discard ~ 2xtau and thin by tau/2 # Adjust the MCMC sampler chains to discard ~ 2xtau and thin by tau/2
# tau = sampler.get_autocorr_time() # tau = sampler.get_autocorr_time()
# print(tau) # print(tau)
result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0] result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0]
fig = plot_result_pdfs(result_dict, xlim=(1, 1.5) if exp == 2 else None) fig = plot_result_pdfs(result_dict, xlim=(1, 1.5) if exp == 2 else None)
fig.savefig(experiment + "result.png") fig.savefig(experiment + "result.pdf")
print(f"posterior mean: {np.mean(rejection_result['Y_accept']):.4e}") print(f"posterior mean: {np.mean(rejection_result['Y_accept']):.4e}")
print(f"posterior uncertainty: {np.std(rejection_result['Y_accept'], ddof=1):.4e}") print(f"posterior uncertainty: {np.std(rejection_result['Y_accept'], ddof=1):.4e}")
......
...@@ -152,7 +152,7 @@ def main(): ...@@ -152,7 +152,7 @@ def main():
fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"]) fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"])
# plt.title(f"Y mean offset: {m_off:.4g}") # plt.title(f"Y mean offset: {m_off:.4g}")
fig.savefig(experiment + "x_samples.png") fig.savefig(experiment + "x_samples.pdf")
result_dict = { result_dict = {
"PRIOR": model["variables"]["Y"]["samples"], "PRIOR": model["variables"]["Y"]["samples"],
"REJECTION": rejection_result["Y_accept"], "REJECTION": rejection_result["Y_accept"],
...@@ -163,14 +163,14 @@ def main(): ...@@ -163,14 +163,14 @@ def main():
sampler = mcmc_reference(model, nwalker=32, nsamples=10000) sampler = mcmc_reference(model, nwalker=32, nsamples=10000)
fig = plot_mcmc_chains(sampler.get_chain(), labels=[key for key, _ in model["variables"].items()]) fig = plot_mcmc_chains(sampler.get_chain(), labels=[key for key, _ in model["variables"].items()])
fig.savefig(experiment + "mcmc.png") fig.savefig(experiment + "mcmc.pdf")
# Adjust the MCMC sampler chains to discard ~ 2xtau and thin by tau/2 # Adjust the MCMC sampler chains to discard ~ 2xtau and thin by tau/2
# tau = sampler.get_autocorr_time() # tau = sampler.get_autocorr_time()
# print(tau) # print(tau)
result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0] result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0]
fig = plot_result_pdfs(result_dict, xlim=(-100, 100)) fig = plot_result_pdfs(result_dict, xlim=(-100, 100))
fig.savefig(experiment + "result.png") fig.savefig(experiment + "result.pdf")
if __name__ == "__main__": if __name__ == "__main__":
......
...@@ -181,7 +181,7 @@ def main(): ...@@ -181,7 +181,7 @@ def main():
fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"], xlim=((8e-6)+1, (13e-6)+1)) fig = plot_x(model["xbar"], model["s"], rejection_result["X_samples"], xlim=((8e-6)+1, (13e-6)+1))
# plt.title(f"Y mean offset: {m_off:.4g}") # plt.title(f"Y mean offset: {m_off:.4g}")
fig.savefig(experiment_run + "x_samples.png") fig.savefig(experiment_run + "x_samples.pdf")
result_dict = { result_dict = {
"PRIOR": model["variables"]["Y"]["samples"], "PRIOR": model["variables"]["Y"]["samples"],
"REJECTION": rejection_result["Y_accept"], "REJECTION": rejection_result["Y_accept"],
...@@ -193,14 +193,14 @@ def main(): ...@@ -193,14 +193,14 @@ def main():
fig = plot_mcmc_chains(sampler.get_chain(), fig = plot_mcmc_chains(sampler.get_chain(),
labels=[key for key, _ in model["variables"].items()]) labels=[key for key, _ in model["variables"].items()])
fig.savefig(experiment_run + "mcmc.png") fig.savefig(experiment_run + "mcmc.pdf")
# Adjust the MCMC sampler chains to discard ~ 2xtau and thin by tau/2 # Adjust the MCMC sampler chains to discard ~ 2xtau and thin by tau/2
# to obtain approximately independent samples. # to obtain approximately independent samples.
# tau = sampler.get_autocorr_time() # tau = sampler.get_autocorr_time()
# print(tau) # print(tau)
result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0] result_dict["MCMC"] = sampler.get_chain(discard=100, flat=True)[:, 0]
fig = plot_result_pdfs(result_dict) fig = plot_result_pdfs(result_dict)
fig.savefig(experiment_run + "result.png") fig.savefig(experiment_run + "result.pdf")
if __name__ == "__main__": if __name__ == "__main__":
......
...@@ -28,14 +28,19 @@ from typing import Any, Union, Optional ...@@ -28,14 +28,19 @@ from typing import Any, Union, Optional
import matplotlib.pyplot as plt # type: ignore import matplotlib.pyplot as plt # type: ignore
import numpy as np import numpy as np
from scipy.stats import norm # type: ignore from scipy.stats import norm # type: ignore
from utility import get_pdf_range, get_hist_data, get_pdf_range_minmax from utility import (get_hist_data, get_pdf_range_minmax)
try:
plt.style.use("code/plot_style.txt")
except IOError:
print("Warning: no plot style is found!")
def plot_x(xbar: float, def plot_x(xbar: float,
std_unc: float, std_unc: float,
x_prior: Union[list, np.ndarray], x_prior: Union[list, np.ndarray],
xlim: Optional[Union[tuple, list]] = None, xlim: Optional[Union[tuple, list]] = None,
ylim: Optional[Union[tuple, list]] = None): ylim: Optional[Union[tuple, list]] = None,
fontsize: Optional[int] = None):
""" """
Plot the prior and data distribution according to given xbar and std_unc, Plot the prior and data distribution according to given xbar and std_unc,
respectively given the samples in x_prior respectively given the samples in x_prior
...@@ -52,10 +57,11 @@ def plot_x(xbar: float, ...@@ -52,10 +57,11 @@ def plot_x(xbar: float,
Returns: Returns:
[type] -- [description] [type] -- [description]
""" """
fig = plt.figure(figsize=(6, 4)) # pylint: disable=too-many-arguments
fig, ax1 = plt.subplots()
# plot X # plot X
hist_grid = get_pdf_range(x_prior, k=2.5, num=1000) # hist_grid = get_pdf_range(x_prior, k=2.5, num=1000)
x_grid, p_s1 = get_hist_data(x_prior, hist_grid) # x_grid, p_s1 = get_hist_data(x_prior, hist_grid)
# plot accepted x # plot accepted x
# p_s1, x_grid = get_hist_data(X_rejection, hist_grid) # p_s1, x_grid = get_hist_data(X_rejection, hist_grid)
...@@ -68,25 +74,29 @@ def plot_x(xbar: float, ...@@ -68,25 +74,29 @@ def plot_x(xbar: float,
minmaxgrid = np.linspace(xlim[0], xlim[1], num=1000) minmaxgrid = np.linspace(xlim[0], xlim[1], num=1000)
data_pdf = norm.pdf(minmaxgrid, loc=xbar, scale=std_unc) data_pdf = norm.pdf(minmaxgrid, loc=xbar, scale=std_unc)
plt.plot(minmaxgrid, data_pdf, label=r"Data N($\bar{x}, s^2$)", color="red", linestyle="solid") ax1.plot(minmaxgrid, data_pdf, label=r"Data N($\bar{x}, s^2$)", color="red", linestyle="solid")
x_grid, p_s1 = get_hist_data(x_prior, minmaxgrid) x_grid, p_s1 = get_hist_data(x_prior, minmaxgrid)
plt.plot(x_grid, p_s1, '-g', label="X prior " + r"$\pi(x)$") ax1.plot(x_grid, p_s1, '-g', label="X prior " + r"$\pi(x)$")
# plt.vlines(np.percentile(X_prior, [2.5, 97.5]), 0, np.max(data_pdf), # plt.vlines(np.percentile(X_prior, [2.5, 97.5]), 0, np.max(data_pdf),
# linestyles="dashed", colors="blue", label="Prior 95% CI") # linestyles="dashed", colors="blue", label="Prior 95% CI")
plt.legend(fontsize=12) if fontsize is not None:
plt.xlabel("type A quantity X [a.u.]", fontsize=14) ax1.legend(fontsize=fontsize)
plt.ylabel("Probability density", fontsize=14) ax1.set_xlabel("type A quantity X [a.u.]", fontsize=fontsize)
ax1.set_ylabel("Probability density", fontsize=fontsize)
else:
ax1.legend()
ax1.set_xlabel("type A quantity X [a.u.]")
ax1.set_ylabel("Probability density")
under = [] under = []
for lia, value1 in enumerate(p_s1): for lia, value1 in enumerate(p_s1):
under.append(np.min([value1, data_pdf[lia]])) under.append(np.min([value1, data_pdf[lia]]))
plt.fill_between(x_grid, [0]*len(under), under, color="blue", alpha=0.5) ax1.fill_between(x_grid, [0]*len(under), under, color="blue", alpha=0.5)
if xlim is not None: if xlim is not None:
plt.xlim(xlim) ax1.set_xlim(xlim)
if ylim is not None: if ylim is not None:
plt.ylim(ylim) ax1.set_ylim(ylim)
plt.tight_layout() plt.tight_layout()
return fig return fig
...@@ -113,6 +123,7 @@ def plot_mcmc_chains(samples: np.ndarray, ...@@ -113,6 +123,7 @@ def plot_mcmc_chains(samples: np.ndarray,
loc_ax.yaxis.set_label_coords(-0.1, 0.5) loc_ax.yaxis.set_label_coords(-0.1, 0.5)
axes[-1].set_xlabel("step number") axes[-1].set_xlabel("step number")
plt.tight_layout()
return fig return fig
...@@ -138,7 +149,7 @@ def plot_result_pdfs(sampler_dict: dict, ...@@ -138,7 +149,7 @@ def plot_result_pdfs(sampler_dict: dict,
Any -- Figure for later plotting or saving Any -- Figure for later plotting or saving
""" """
assert "REJECTION" in sampler_dict assert "REJECTION" in sampler_dict
fig = plt.figure(figsize=(6, 4)) fig, ax1 = plt.subplots()
# get grid to plot on from all given densities # get grid to plot on from all given densities
if xlim is not None: if xlim is not None:
...@@ -147,31 +158,31 @@ def plot_result_pdfs(sampler_dict: dict, ...@@ -147,31 +158,31 @@ def plot_result_pdfs(sampler_dict: dict,
hist_grid = get_pdf_range_minmax(sampler_dict, k=5, num=100) hist_grid = get_pdf_range_minmax(sampler_dict, k=5, num=100)
# rejection posterior # rejection posterior
plt.plot(*get_hist_data(sampler_dict["REJECTION"], hist_grid), '-b', label="Rejection") ax1.plot(*get_hist_data(sampler_dict["REJECTION"], hist_grid), '-b', label="Rejection")
# prior # prior
if "PRIOR" in sampler_dict: if "PRIOR" in sampler_dict:
plt.plot(*get_hist_data(sampler_dict["PRIOR"], hist_grid), ax1.plot(*get_hist_data(sampler_dict["PRIOR"], hist_grid),
linestyle="dashed", color="black", label="Prior") linestyle="dashed", color="black", label="Prior")
# S1 posterior # S1 posterior
if "GUMS1" in sampler_dict: if "GUMS1" in sampler_dict:
plt.plot(*get_hist_data(sampler_dict["GUMS1"], hist_grid), '-g', label="GUM-S1") ax1.plot(*get_hist_data(sampler_dict["GUMS1"], hist_grid), '-g', label="GUM-S1")
# MCMC posterior # MCMC posterior
if "MCMC" in sampler_dict: if "MCMC" in sampler_dict:
plt.hist(sampler_dict["MCMC"], ax1.hist(sampler_dict["MCMC"],
bins=100, label="MCMC", density=True, alpha=0.5, color="blue", bins=100, label="MCMC", density=True, alpha=0.5, color="blue",
edgecolor="blue") edgecolor="blue")
# GUM posterior # GUM posterior
if "GUM" in sampler_dict: if "GUM" in sampler_dict:
plt.plot(hist_grid, norm.pdf(hist_grid, loc=sampler_dict["GUM"]["mean"], ax1.plot(hist_grid, norm.pdf(hist_grid, loc=sampler_dict["GUM"]["mean"],
scale=sampler_dict["GUM"]["uncertainty"]), '-r', label="GUM") scale=sampler_dict["GUM"]["uncertainty"]), '-r', label="GUM")
plt.legend(fontsize=12) ax1.legend(fontsize=12)
# plt.title(f"Y mean offset: {m_off:.4g}, Y noise faktor: {s_off:.2g}", fontsize=14) # plt.title(f"Y mean offset: {m_off:.4g}, Y noise faktor: {s_off:.2g}", fontsize=14)
if xlim is not None: if xlim is not None:
plt.xlim(xlim) plt.xlim(xlim)
plt.xlabel("Y", fontsize=14) ax1.set_xlabel("measurand Y [a.u.]")
plt.ylabel("Probability density", fontsize=14) ax1.set_ylabel("Probability density")
plt.tight_layout()
return fig return fig
...@@ -28,6 +28,10 @@ import os ...@@ -28,6 +28,10 @@ import os
import numpy as np import numpy as np
import matplotlib.pyplot as plt # type: ignore import matplotlib.pyplot as plt # type: ignore
from scipy.stats import norm # type: ignore from scipy.stats import norm # type: ignore
try:
plt.style.use("code/plot_style.txt")
except IOError:
print("Warning: no plot style is found!")
def main(n_samples: int): def main(n_samples: int):
...@@ -59,18 +63,19 @@ def main(n_samples: int): ...@@ -59,18 +63,19 @@ def main(n_samples: int):
rejectedx.append(q_samples[lia]) rejectedx.append(q_samples[lia])
rejectedy.append(uniform_sample) rejectedy.append(uniform_sample)
fig = plt.figure() fig, ax1 = plt.subplots()
plt.plot(x_grid, q_den.pdf(x_grid), label="prior") ax1.plot(x_grid, q_den.pdf(x_grid), label="prior")
plt.plot(x_grid, p_den.pdf(x_grid), label="posterior") ax1.plot(x_grid, p_den.pdf(x_grid), label="posterior")
plt.plot(x_grid, supremum*q_den.pdf(x_grid), linestyle="dashed", color="black", label="envelope") ax1.plot(x_grid, supremum*q_den.pdf(x_grid), linestyle="dashed", color="black", label="envelope")
plt.scatter(acceptedx, acceptedy, marker="o", color="green", s=2, alpha=0.5, edgecolors="green") ax1.scatter(acceptedx, acceptedy, marker="o", color="green", s=2, alpha=0.5, edgecolors="green")
plt.scatter(rejectedx, rejectedy, marker="o", color="red", s=2, alpha=0.5, edgecolors="red") ax1.scatter(rejectedx, rejectedy, marker="o", color="red", s=2, alpha=0.5, edgecolors="red")
plt.ylabel("Probability density function") ax1.set_ylabel("Probability density function")
plt.xlabel("Y") ax1.set_xlabel("measurand Y [a.u.]")
plt.xlim(-10, 11) ax1.set_xlim(-10, 11)
plt.title(f"Number of samples: N={n_samples}") # ax.set_title(f"Number of samples: N={n_samples}")
plt.legend() ax1.legend()
fig.savefig(f"data/sampler_visualization/N_{n_samples}.pdf", dpi=300) plt.tight_layout()
fig.savefig(f"data/sampler_visualization/N_{n_samples}.pdf")
if __name__ == "__main__": if __name__ == "__main__":
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment