From b2f6c58490cd3517d1a813e1b51b4e2aafc945a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Patrik=20M=C3=BCller?= Date: Thu, 2 Jun 2022 08:09:02 +0200 Subject: Add relative error plots with shift --- buch/papers/laguerre/scripts/gamma_approx.py | 97 ++++++++++++++++++++++++---- 1 file changed, 84 insertions(+), 13 deletions(-) (limited to 'buch/papers/laguerre/scripts/gamma_approx.py') diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py index 90843b1..9c8f3ee 100644 --- a/buch/papers/laguerre/scripts/gamma_approx.py +++ b/buch/papers/laguerre/scripts/gamma_approx.py @@ -37,11 +37,42 @@ def find_shift(z, target): return zs, factor -def laguerre_gamma_shift(z, x=None, w=None, n=8, target=11): +def find_optimal_shift(z, n): + mhat = 1.34093 * n + 0.854093 + steps = int(np.ceil(mhat - np.real(z))) - 1 + return steps + + +def get_shifting_factor(z, steps): + if steps > 0: + factor = 1 / pochhammer(z, steps) + elif steps < 0: + factor = pochhammer(z + steps, -steps) + return factor + + +def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11): x, w = _prep_zeros_and_weights(x, w, n) + n = len(x) z += 0j z_shifted, correction_factor = find_shift(z, target) + + res = np.sum(x ** (z_shifted - 1) * w) + res *= correction_factor + res = drop_imag(res) + return res + + +def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8): + x, w = _prep_zeros_and_weights(x, w, n) + n = len(x) + + z += 0j + opt_shift = find_optimal_shift(z, n) + correction_factor = get_shifting_factor(z, opt_shift) + z_shifted = z + opt_shift + res = np.sum(x ** (z_shifted - 1) * w) res *= correction_factor res = drop_imag(res) @@ -72,8 +103,10 @@ def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs): f = laguerre_gamma_simple elif func == "mirror": f = laguerre_gamma_mirror + elif func == "optimal_shifted": + f = laguerre_gamma_opt_shifted else: - f = laguerre_gamma_shift + f = laguerre_gamma_shifted return np.array([f(zi, x, w, n, **kwargs) for zi in z]) @@ -81,11 +114,10 @@ def calc_rel_error(x, y): return (y - x) / x -ns = np.arange(2, 12, 2) - # Simple / naive xmin = -5 xmax = 30 +ns = np.arange(2, 12, 2) ylim = np.array([-11, 6]) x = np.linspace(xmin + EPSILON, xmax - EPSILON, 400) gamma = scipy.special.gamma(x) @@ -104,7 +136,7 @@ ax.set_xlabel(r"$z$") ax.set_ylabel("Relativer Fehler") ax.legend(ncol=3, fontsize="small") ax.grid(1, "both") -fig.savefig(f"{img_path}/rel_error_simple.pgf") +# fig.savefig(f"{img_path}/rel_error_simple.pgf") # Mirrored @@ -130,7 +162,7 @@ ax2.set_xlabel(r"$z$") ax2.set_ylabel("Relativer Fehler") ax2.legend(ncol=1, loc="upper left", fontsize="small") ax2.grid(1, "both") -fig2.savefig(f"{img_path}/rel_error_mirror.pgf") +# fig2.savefig(f"{img_path}/rel_error_mirror.pgf") # Move to target @@ -163,12 +195,14 @@ v = ax3.imshow(bests, cmap="inferno", aspect="auto", interpolation="nearest") plt.colorbar(v, ax=ax3, label=r"$m$") ticks = np.arange(0, N + 1, N // 5) ax3.set_xlim(0, 1) -ax3.set_xticks(ticks, [f"{v:.2f}" for v in ticks / N]) +ax3.set_xticks(ticks) +ax3.set_xticklabels([f"{v:.2f}" for v in ticks / N]) ax3.set_xticks(np.arange(0, N + 1, N // 20), minor=True) -ax3.set_yticks(np.arange(len(ns)), ns) +ax3.set_yticks(np.arange(len(ns))) +ax3.set_yticklabels(ns) ax3.set_xlabel(r"$z$") ax3.set_ylabel(r"$n$") -fig3.savefig(f"{img_path}/targets.pdf") +# fig3.savefig(f"{img_path}/targets.pdf") targets = np.mean(bests, -1) intercept, bias = np.polyfit(ns, targets, 1) @@ -178,9 +212,7 @@ fig4, axs4 = plt.subplots( xl = np.array([ns[0] - 0.5, ns[-1] + 0.5]) axs4[0].plot(xl, intercept * xl + bias, label=r"$\hat{m}$") axs4[0].plot(ns, targets, "x", label=r"$\bar{m}$") -axs4[1].plot( - ns, ((intercept * ns + bias) - targets), "-x", label=r"$\hat{m} - \bar{m}$" -) +axs4[1].plot(ns, ((intercept * ns + bias) - targets), "-x", label=r"$\hat{m} - \bar{m}$") axs4[0].set_xlim(*xl) # axs4[0].set_title("Schätzung von Mittelwert") # axs4[1].set_title("Fehler") @@ -188,10 +220,49 @@ axs4[-1].set_xlabel(r"$z$") for ax in axs4: ax.grid(1) ax.legend() -fig4.savefig(f"{img_path}/schaetzung.pgf") +# fig4.savefig(f"{img_path}/schaetzung.pgf") print(f"Intercept={intercept:.6g}, Bias={bias:.6g}") predicts = np.ceil(intercept * ns[:, None] + bias - x) print(f"Error: {int(np.sum(np.abs(bests-predicts)))}") +# Comparison relative error between methods +N = 200 +step = 1 / (N - 1) +x = np.linspace(step, 1 - step, N + 1) +gamma = scipy.special.gamma(x)[:, None] +n = 8 +targets = np.arange(10, 14) +gamma = scipy.special.gamma(x) +fig5, ax5 = plt.subplots(num=1, clear=True, constrained_layout=True) +for target in targets: + gamma_lag = eval_laguerre_gamma(x, target=target, n=n, func="shifted") + rel_error = np.abs(calc_rel_error(gamma, gamma_lag)) + ax5.semilogy(x, rel_error, label=f"$m={target}$") +gamma_lgo = eval_laguerre_gamma(x, n=n, func="optimal_shifted") +rel_error = np.abs(calc_rel_error(gamma, gamma_lgo)) +ax5.semilogy(x, rel_error, label="$m^*$") +ax5.set_xlim(x[0], x[-1]) +ax5.set_ylim(5e-9, 5e-8) +ax5.set_xlabel(r"$z$") +ax5.grid(1, "both") +ax5.legend() +fig5.savefig(f"{img_path}/rel_error_shifted.pgf") + +N = 200 +x = np.linspace(-5+ EPSILON, 5-EPSILON, N) +gamma = scipy.special.gamma(x)[:, None] +n = 8 +gamma = scipy.special.gamma(x) +fig6, ax6 = plt.subplots(num=1, clear=True, constrained_layout=True) +gamma_lgo = eval_laguerre_gamma(x, n=n, func="optimal_shifted") +rel_error = np.abs(calc_rel_error(gamma, gamma_lgo)) +ax6.semilogy(x, rel_error, label="$m^*$") +ax6.set_xlim(x[0], x[-1]) +ax6.set_ylim(5e-9, 5e-8) +ax6.set_xlabel(r"$z$") +ax6.grid(1, "both") +ax6.legend() +fig6.savefig(f"{img_path}/rel_error_range.pgf") + # plt.show() -- cgit v1.2.1