commit b78756847872f6eb528b4ba893d807c7a2dfbfa0 parent 0b8b6816f2827367014542a4024b74df7d5cabfe Author: miksa <milutin@popovic.xyz> Date: Sun, 11 Apr 2021 14:31:25 +0200 chi^2 dof test implemented with plots Diffstat:
24 files changed, 95 insertions(+), 62 deletions(-)
diff --git a/sesh1/prog/t0-method/__pycache__/method.cpython-39.pyc b/sesh1/prog/t0-method/__pycache__/method.cpython-39.pyc Binary files differ. diff --git a/sesh1/prog/t0-method/__pycache__/model.cpython-39.pyc b/sesh1/prog/t0-method/__pycache__/model.cpython-39.pyc Binary files differ. diff --git a/sesh1/prog/t0-method/method.py b/sesh1/prog/t0-method/method.py @@ -5,34 +5,35 @@ import sympy as sp global x; x = sp.symbols("x") -def design_matrix(model, var): - A = sp.diff(model(x, *var), var.T) - A_f = sp.lambdify((x, *var), A, modules=['numpy', {'Heaviside': np.heaviside}]) - return A_f +def t0_fit(model, var_str, x_data, y_data,\ + p0, cov_stat, cov_relsyst, way=0, iterations=100, alpha=0.1): + var = sp.Matrix(sp.symbols(var_str)) -def f2_lambda(model, var): - f2_f = sp.lambdify((x, *var), model(x, *var), modules=['numpy', {'Heaviside': np.heaviside}]) - return f2_f + jacobi = sp.diff(model(x, *var), var.T) + A_func = sp.lambdify((x, *var), jacobi, modules=['numpy', {'Heaviside': np.heaviside}]) -def t0_fit(model, var_str, x_data, y_data,\ - p0, cov_stat, cov_relsyst, iterations=100, alpha=0.1): + f2_func = sp.lambdify((x, *var), model(x, *var), modules=['numpy', {'Heaviside': np.heaviside}]) - var = sp.Matrix(sp.symbols(var_str)) - A_func = design_matrix(model, var) - f2_func = f2_lambda(model, var) + chi_sq = [] for i in range(iterations): A = A_func(x_data, *p0)[0].T f2 = f2_func(x_data, *p0) - P = np.linalg.inv(cov_relsyst @ np.outer(f2, f2) + cov_stat) + if way == 0: # consider D'Agostini bias + P = np.linalg.inv(cov_relsyst @ np.outer(f2, f2) + cov_stat) + else: # don't consider D'Agostini bias + P = np.linalg.inv(cov_relsyst * (p0@p0) + cov_stat) delta_p = np.linalg.inv(A.T @ P @ A) @ A.T @ P @ (y_data - f2) p0 = p0 + alpha*delta_p + + chi_sq.append((y_data - f2).T @ P @ (y_data - f2)) + print(f"{i} Iterations completed") print('Done\n') dp = np.sqrt(np.diag(np.linalg.inv(A.T@P@A))) - return p0, dp + return p0, dp, np.array(chi_sq) diff --git a/sesh1/prog/t0-method/model.py b/sesh1/prog/t0-method/model.py @@ -4,7 +4,6 @@ import numpy as np import sympy as sp import matplotlib.pyplot as plt - global m_p; m_p = 0.13957 sig_p = lambda x: sp.sqrt(1 - 4*m_p**2/x) @@ -16,37 +15,3 @@ def model(s, m_q, g_q, m_w, g_w, e_w, a, b, c): part3 = (1 + a*s + b*s**2 + c*s**3)**2 return part1 * part2 * part3 - -def my_plot(name, x_data, y_data, p, dp, sigma): - - for i in range(len(p)): - print(f'{p[i]}\t\\pm {dp[i]}') - - x_model = np.linspace(x_data[0], x_data[-1], 500) - - x = sp.symbols('x') - la_mod = model(x, *p) - la_mod = sp.lambdify(x, la_mod, modules=['numpy', {'Heaviside': np.heaviside}]) - y_model = la_mod(x_model) - - plt.figure(figsize=[10, 7]) - plt.errorbar(x_data, y_data, yerr=sigma, c='black', fmt='.k', label=f'{name}') - plt.plot(x_model, y_model, c='red', label='Fit') - - p, dp = np.round(p, 5), np.round(dp, 3) - plt.annotate(r'$M_{\rho} = $' + f'({p[0]}' + r'$\pm$' + f'{dp[0]}) GeV', (0.2, 40)) - plt.annotate(r'$\Gamma_{\rho} = $' + f'({p[1]}' + r'$\pm$' + f'{dp[1]}) GeV', (0.2, 38)) - plt.annotate(r'$M_{\omega} = $' + f'({p[2]}' + r'$\pm$' + f'{dp[2]}) GeV', (0.2, 36)) - plt.annotate(r'$\Gamma_{\omega} = $' + f'({p[3]}' + r'$\pm$' + f'{dp[3]}) GeV', (0.2, 34)) - plt.annotate(r'$\epsilon_{\omega} = $' + f'({p[4]}' + r'$\pm$' + f'{dp[4]}) GeV', (0.2, 32)) - plt.annotate(r'$a = $' + f'({p[5]}' + r'$\pm$' + f'{dp[5]}) GeV', (0.2, 30)) - plt.annotate(r'$b = $' + f'({p[6]}' + r'$\pm$' + f'{dp[6]}) GeV', (0.2, 28)) - plt.annotate(r'$c = $' + f'({p[7]}' + r'$\pm$' + f'{dp[7]}) GeV', (0.2, 26)) - - plt.title(f't0-Singlefit of {name}') - plt.legend(loc='best') - plt.ylabel(r'$|F_{\pi}^V(s)|^2|$') - plt.xlabel(r'$\sqrt{s}$[GeV]') - - plt.savefig(f'./plots/{name}.png') - plt.close() diff --git a/sesh1/prog/t0-method/multiple_fit.py b/sesh1/prog/t0-method/multiple_fit.py @@ -3,6 +3,7 @@ import numpy as np import sympy as sp from scipy.linalg import block_diag from matplotlib import pyplot as plt +from scipy.special import gammaincc from method import * from model import * @@ -52,7 +53,7 @@ def main(): var_str = "m_q g_q m_w g_w e_w a b c" # Fit - p, dp = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ + p, dp, chi_sq = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ cov_relsyst) @@ -91,5 +92,18 @@ def main(): plt.savefig('./plots/all-fit.png') plt.close() + dof = len(x_data) - len(p) + p_value = 1 - gammaincc(chi_sq/2, dof/2) + + plt.figure(figsize=[10, 7]) + plt.plot(chi_sq/dof, p_value, label=r'$\chi^2_{min}/dof = $' + f'{round(chi_sq[-1]/dof, 3)}') + plt.title(r'$\chi^2$'+'Multifit') + plt.xlabel(r'$\chi^2/dof$') + plt.ylabel('p-value') + plt.legend(loc='best') + + plt.savefig(f'./plots/all-chisq.png') + plt.close() + if __name__ == "__main__": main() diff --git a/sesh1/prog/t0-method/my_plot.py b/sesh1/prog/t0-method/my_plot.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3.9 + +import numpy as np +import sympy as sp +import matplotlib.pyplot as plt +from scipy.special import gammaincc + +from model import * + +def my_plot(name, x_data, y_data, p, dp, chi_sq, sigma): + + for i in range(len(p)): + print(f'{p[i]}\t\\pm {dp[i]}') + + x_model = np.linspace(x_data[0], x_data[-1], 500) + + x = sp.symbols('x') + la_mod = model(x, *p) + la_mod = sp.lambdify(x, la_mod, modules=['numpy', {'Heaviside': np.heaviside}]) + y_model = la_mod(x_model) + + plt.figure(figsize=[10, 7]) + plt.errorbar(x_data, y_data, yerr=sigma, c='black', fmt='.k', label=f'{name}') + plt.plot(x_model, y_model, c='red', label='Fit') + + p, dp = np.round(p, 5), np.round(dp, 3) + plt.annotate(r'$M_{\rho} = $' + f'({p[0]}' + r'$\pm$' + f'{dp[0]}) GeV', (0.2, 40)) + plt.annotate(r'$\Gamma_{\rho} = $' + f'({p[1]}' + r'$\pm$' + f'{dp[1]}) GeV', (0.2, 38)) + plt.annotate(r'$M_{\omega} = $' + f'({p[2]}' + r'$\pm$' + f'{dp[2]}) GeV', (0.2, 36)) + plt.annotate(r'$\Gamma_{\omega} = $' + f'({p[3]}' + r'$\pm$' + f'{dp[3]}) GeV', (0.2, 34)) + plt.annotate(r'$\epsilon_{\omega} = $' + f'({p[4]}' + r'$\pm$' + f'{dp[4]}) GeV', (0.2, 32)) + plt.annotate(r'$a = $' + f'({p[5]}' + r'$\pm$' + f'{dp[5]}) GeV', (0.2, 30)) + plt.annotate(r'$b = $' + f'({p[6]}' + r'$\pm$' + f'{dp[6]}) GeV', (0.2, 28)) + plt.annotate(r'$c = $' + f'({p[7]}' + r'$\pm$' + f'{dp[7]}) GeV', (0.2, 26)) + + plt.title(f't0-Singlefit of {name}') + plt.legend(loc='best') + plt.ylabel(r'$|F_{\pi}^V(s)|^2|$') + plt.xlabel(r'$\sqrt{s}$[GeV]') + + plt.savefig(f'./plots/{name}.png') + plt.close() + + + dof = len(x_data) - len(p) + p_value = 1 - gammaincc(chi_sq/2, dof/2) + + plt.figure(figsize=[10, 7]) + plt.plot(chi_sq/dof, p_value, label=r'$\chi^2_{min}/dof = $' + f'{round(chi_sq[-1]/dof, 3)}') + plt.title(r'$\chi^2$'+f'-{name}') + plt.xlabel(r'$\chi^2/dof$') + plt.ylabel('p-value') + plt.legend(loc='best') + + plt.savefig(f'./plots/{name}-chisq.png') + plt.close() diff --git a/sesh1/prog/t0-method/plots/BABAR-KLOE-chisq.png b/sesh1/prog/t0-method/plots/BABAR-KLOE-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/BABAR-KLOE-fit.png b/sesh1/prog/t0-method/plots/BABAR-KLOE-fit.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/BABAR-chisq.png b/sesh1/prog/t0-method/plots/BABAR-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/CMD2-BABAR-chisq.png b/sesh1/prog/t0-method/plots/CMD2-BABAR-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/CMD2-BABAR-fit.png b/sesh1/prog/t0-method/plots/CMD2-BABAR-fit.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/CMD2-KLOE-chisq.png b/sesh1/prog/t0-method/plots/CMD2-KLOE-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/CMD2-KLOE-fit.png b/sesh1/prog/t0-method/plots/CMD2-KLOE-fit.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/CMD2-chisq.png b/sesh1/prog/t0-method/plots/CMD2-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/KLOE-chisq.png b/sesh1/prog/t0-method/plots/KLOE-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/SND-BABAR-chisq.png b/sesh1/prog/t0-method/plots/SND-BABAR-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/SND-BABAR-fit.png b/sesh1/prog/t0-method/plots/SND-BABAR-fit.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/SND-CMD2-chisq.png b/sesh1/prog/t0-method/plots/SND-CMD2-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/SND-CMD2-fit.png b/sesh1/prog/t0-method/plots/SND-CMD2-fit.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/SND-KLOE-chisq.png b/sesh1/prog/t0-method/plots/SND-KLOE-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/SND-KLOE-fit.png b/sesh1/prog/t0-method/plots/SND-KLOE-fit.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/SND-chisq.png b/sesh1/prog/t0-method/plots/SND-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/plots/all-chisq.png b/sesh1/prog/t0-method/plots/all-chisq.png Binary files differ. diff --git a/sesh1/prog/t0-method/single_fit.py b/sesh1/prog/t0-method/single_fit.py @@ -1,9 +1,6 @@ #!/usr/bin/env python3.9 -import numpy as np -import sympy as sp -from matplotlib import pyplot as plt - +from my_plot import * from method import * from model import * @@ -18,11 +15,11 @@ def SND(): cov_relsyst = (np.eye(len(x_data)) * np.array([0.032 if i<=2 else 0.013 for i in range(len(x_data))]))**2 var_str = "m_q g_q m_w g_w e_w a b c" - p, dp = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ + p, dp, chi_sq = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ cov_relsyst) sigma = np.diag(cov_stat) - my_plot('SND', x_data, y_data, p, dp, sigma) + my_plot('SND', x_data, y_data, p, dp, chi_sq, sigma) def CMD2(): data = np.loadtxt('../data/CMD2-VFF.txt') @@ -42,11 +39,11 @@ def CMD2(): var_str = "m_q g_q m_w g_w e_w a b c" - p, dp = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ + p, dp, chi_sq = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ cov_relsyst) sigma = np.sqrt(np.diag(cov_stat)) - my_plot('CMD2', x_data, y_data, p, dp, sigma) + my_plot('CMD2', x_data, y_data, p, dp, chi_sq, sigma) def KLOE(): data = np.loadtxt('../data/KLOE-VFF.txt') @@ -56,11 +53,11 @@ def KLOE(): cov_stat = np.loadtxt('../data/KLOE-StatCov.txt') var_str = "m_q g_q m_w g_w e_w a b c" - p, dp = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ + p, dp, chi_sq = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ cov_relsyst) sigma = np.sqrt(np.diag(cov_stat)) - my_plot('KLOE', x_data, y_data, p, dp, sigma) + my_plot('KLOE', x_data, y_data, p, dp, chi_sq, sigma) def BABAR(): data = np.loadtxt('../data/BABAR-VFF.txt') @@ -70,11 +67,11 @@ def BABAR(): cov_stat = np.loadtxt('../data/BABAR-StatCov.txt') var_str = "m_q g_q m_w g_w e_w a b c" - p, dp = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ + p, dp, chi_sq = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\ cov_relsyst) sigma = np.sqrt(np.diag(cov_stat)) - my_plot('BABAR', x_data, y_data, p, dp, sigma) + my_plot('BABAR', x_data, y_data, p, dp, chi_sq, sigma) def main(): SND()
