Skip to content
Snippets Groups Projects
Commit 5aeb44a5 authored by Hartmut Stadie's avatar Hartmut Stadie
Browse files

use only notebook, fixes

parent 67874142
No related branches found
No related tags found
No related merge requests found
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "3885a60f-ddf0-4c7b-a587-6e7e037ec9b5",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.stats as stats\n",
"import matplotlib.pyplot as plt \n",
"\n",
"def f(x):\n",
" return 2*x + 1\n",
"\n",
"n = 10\n",
"xs = np.linspace(0,4,n)\n",
"sigma_y=0.4\n",
"ys = stats.multivariate_normal.rvs(f(xs), np.eye(n)*sigma_y**2, 1, random_state=42)\n",
" \n",
"x_axis = np.linspace(0,4,100)\n",
"plt.errorbar(xs,ys,yerr=sigma_y,fmt=\".\")\n",
"plt.plot(x_axis, f(x_axis),'--')\n",
"plt.xlabel(\"x\")\n",
"plt.ylabel(\"y\")\n",
"plt.savefig(\"line.png\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9f78c037-4cca-448c-829c-dac9dea11d3c",
"metadata": {},
"outputs": [],
"source": [
"def chi2(x, y, sy, a, m):\n",
" my = m * x + a\n",
" r = (y - my)/sy\n",
" return np.sum(r**2)\n",
"\n",
"chi2(xs, ys, sigma_y, 1, 2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "76181537-7b07-4329-9615-0947d4b9dd5a",
"metadata": {},
"outputs": [],
"source": [
"from matplotlib import cm\n",
"a_mesh, m_mesh = np.meshgrid(np.linspace(-3,5,100), np.linspace(-3,5,100))\n",
"chi2_vect = np.vectorize(chi2, excluded=[0, 1])\n",
"m_axis = np.linspace(-3,5,100)\n",
"plt.plot(m_axis, chi2_vect(xs, ys, sigma_y, 1, m_axis))\n",
"plt.xlabel(\"m\")\n",
"plt.show()\n",
"plt.plot(m_axis, chi2_vect(xs, ys, sigma_y, m_axis, 2))\n",
"plt.xlabel(\"a\")\n",
"plt.show()\n",
"\n",
"plt.contourf(a_mesh, m_mesh, chi2_vect(xs, ys, sigma_y, a_mesh, m_mesh), vmin=0, vmax=5000, levels=200, cmap=cm.coolwarm)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "34aaaecd-6d6a-452d-8602-d8b7b9d34b72",
"metadata": {},
"source": [
"## Analytisch\n",
"$\\hat m = \\frac{\\langle xy \\rangle - \\langle y \\rangle\\langle x \\rangle}{\\langle x^2 \\rangle - \\langle x \\rangle^2}$\n",
"\n",
"$\\hat a = \\frac{ \\langle y \\rangle \\langle x^2 \\rangle - \\langle x \\rangle \\langle xy \\rangle}{ \\langle x^2 \\rangle - \\langle x \\rangle^2}$\n",
"\n",
"\n",
"$V(\\hat m) = \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{1}{\\langle x^2 \\rangle - \\langle x \\rangle^2}$\n",
"\n",
"\n",
"$V(\\hat a) = \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x^2 \\rangle}{\\langle x^2 \\rangle - \\langle x \\rangle^2}$\n",
"\n",
"\n",
"$\\text{cov}(\\hat m, \\hat a) = - \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{<x>}{\\langle x^2 \\rangle - \\langle x \\rangle^2}$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d3b470d1-1b48-46dc-ae1c-d7a44be51e8c",
"metadata": {},
"outputs": [],
"source": [
"C = np.eye(n)*sigma_y*sigma_y\n",
"W = np.linalg.inv(C)\n",
"J = np.ones((len(ys),1))\n",
"sum_xy = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*ys)))[0]\n",
"sum_y = (np.linalg.inv(J.T@W@J)@(J.T@W@(ys)))[0]\n",
"sum_x = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs)))[0]\n",
"sum_x2 = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*xs)))[0]\n",
"print(sum_xy, sum_y, sum_x, sum_x2)\n",
"mhat = (sum_xy - sum_x*sum_y)/(sum_x2-sum_x*sum_x)\n",
"ahat = (sum_y*sum_x2 - sum_x*sum_xy)/(sum_x2-sum_x*sum_x)\n",
"print(mhat, ahat)\n",
"sum_siginv = sigma_y*sigma_y/len(ys)\n",
"V_am = np.array([[sum_siginv/(sum_x2-sum_x*sum_x), -sum_siginv*sum_x/(sum_x2-sum_x*sum_x)],\n",
" [-sum_siginv*sum_x/(sum_x2-sum_x*sum_x), sum_siginv*sum_x2/(sum_x2-sum_x*sum_x)]])\n",
"print(V_am)\n",
"print(np.sqrt(V_am[0,0]), np.sqrt(V_am[1,1]), V_am[1,0]/(np.sqrt(V_am[0,0])*np.sqrt(V_am[1,1])))"
]
},
{
"cell_type": "markdown",
"id": "2ba98419-a39d-48e5-8ed1-47fdcab7d237",
"metadata": {},
"source": [
"# Kovarianz direkt aus C\n",
"$\\frac{d\\hat m}{y_i} = \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{x_i - \\langle x \\rangle}{\\sigma_i^2(\\langle x^2 \\rangle - \\langle x \\rangle^2)}$\n",
"\n",
"$\\frac{d\\hat a}{y_i} = \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x^2 \\rangle - \\langle x \\rangle x_i}{\\sigma_i^2(\\langle x^2 \\rangle - \\langle x \\rangle^2)}$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bf423c2c-0c8f-429e-ab36-4c2d074d3575",
"metadata": {},
"outputs": [],
"source": [
"dmdy = sum_siginv / (sum_x2-sum_x*sum_x)*(xs-sum_x)/sigma_y**2\n",
"dady = sum_siginv / (sum_x2-sum_x*sum_x)*(sum_x2-xs*sum_x)/sigma_y**2\n",
"A=np.array([dmdy,dady]).T\n",
"print(A)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "58615ab8-b1a8-433e-9380-666185592201",
"metadata": {},
"outputs": [],
"source": [
"C_am = A.T@C@A\n",
"print(C_am)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fe255a0f-c44d-4fc0-ab4b-e4bcb7c16f70",
"metadata": {},
"outputs": [],
"source": [
"x_axis = np.linspace(0,4,100)\n",
"plt.errorbar(xs,ys,yerr=sigma_y,fmt=\".\")\n",
"plt.plot(x_axis, f(x_axis),'--')\n",
"plt.plot(x_axis, x_axis*mhat + ahat,'-')\n",
"plt.xlabel(\"x\")\n",
"plt.ylabel(\"y\")\n",
"plt.savefig(\"line.png\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "b2b90d30-4e37-4ecc-866d-a45b2243d8bc",
"metadata": {},
"source": [
"# Fehler auf Ausgleichsgerade\n",
"\n",
"$y = \\hat m x +\\hat a$ \n",
"\n",
"$\\frac{dy}{dm} = x$ und $\\frac{dy}{da} = 1$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2d61e62e-77c1-41ca-a0e9-2f58bb44de08",
"metadata": {},
"outputs": [],
"source": [
"def err(x, V):\n",
" A = np.array([x, 1])\n",
" return np.sqrt(A.T@V@A)\n",
"\n",
"err(0, V_am)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c2f71b37-b37c-4802-a356-e57d9c8df0de",
"metadata": {},
"outputs": [],
"source": [
"x_axis = np.linspace(0,4,100)\n",
"plt.errorbar(xs,ys,yerr=sigma_y,fmt=\".\")\n",
"plt.plot(x_axis, f(x_axis),'--')\n",
"plt.plot(x_axis, x_axis*mhat + ahat,'-g')\n",
"plt.plot(x_axis, x_axis*mhat + ahat - np.vectorize(err,excluded=[1])(x_axis, V_am),'--g')\n",
"plt.plot(x_axis, x_axis*mhat + ahat + np.vectorize(err,excluded=[1])(x_axis, V_am),'--g')\n",
"\n",
"plt.xlabel(\"x\")\n",
"plt.ylabel(\"y\")\n",
"plt.savefig(\"line.png\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f76530a7-f385-4e58-a504-1f7e75f77129",
"metadata": {},
"outputs": [],
"source": [
"plt.plot(x_axis, np.vectorize(err, excluded=[1])(x_axis, V_am),'-')"
]
},
{
"cell_type": "markdown",
"id": "952ccb17-cf14-4561-ad21-acc50b72cf73",
"metadata": {},
"source": [
"$\\chi^2 = (\\sum_i \\frac{1}{\\sigma_i^2}) V(y) ( 1- \\rho^2_{xy})$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "348c5655-6356-4e5b-8f88-f4cbbd4739f7",
"metadata": {},
"outputs": [],
"source": [
"print(len(ys)/(sigma_y*sigma_y)*np.var(ys)*(1-np.corrcoef(xs,ys)[0,1]**2))"
]
},
{
"cell_type": "markdown",
"id": "e6df3805-e62d-4cba-ba9f-1f0c2f187d55",
"metadata": {},
"source": [
"$\\chi^2 =\\sum \\frac{(y -\\hat m x - \\hat a)^2}{\\sigma_y^2}$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ec6c4161-50dd-475b-8f0f-d02eb9f18cfe",
"metadata": {},
"outputs": [],
"source": [
"np.sum(((ys - (xs * mhat + ahat))/sigma_y)**2)"
]
},
{
"cell_type": "markdown",
"id": "bd398a45-9d44-4a28-8cec-592a9b38ec99",
"metadata": {},
"source": [
"# $\\chi^2$ für Ausgleichsgeraden"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5e56c09a-f030-4dad-a095-6b28614d36b2",
"metadata": {},
"outputs": [],
"source": [
"ys_set = stats.multivariate_normal.rvs(f(xs), np.eye(n)*sigma_y**2, 1000, random_state=42)\n",
"chi2s = np.zeros(len(ys_set))\n",
"for i in range(len(ys_set)):\n",
" sum_xy = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*ys_set[i])))[0]\n",
" sum_y = (np.linalg.inv(J.T@W@J)@(J.T@W@(ys_set[i])))[0]\n",
" sum_x = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs)))[0]\n",
" sum_x2 = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*xs)))[0]\n",
" mhat = (sum_xy - sum_x*sum_y)/(sum_x2-sum_x*sum_x)\n",
" ahat = (sum_y*sum_x2 - sum_x*sum_xy)/(sum_x2-sum_x*sum_x)\n",
" sum_siginv = sigma_y*sigma_y/len(ys)\n",
" V_am = np.array([[sum_siginv/(sum_x2-sum_x*sum_x), -sum_siginv*sum_x/(sum_x2-sum_x*sum_x)],\n",
" [-sum_siginv*sum_x/(sum_x2-sum_x*sum_x), sum_siginv*sum_x2/(sum_x2-sum_x*sum_x)]])\n",
" chi2s[i] = np.sum(((ys_set[i] - (xs * mhat + ahat))/sigma_y)**2)\n",
" #print(mhat, ahat, ys_set[i], xs * mhat + ahat, (ys_set[i] - (xs * mhat + ahat))/sigma_y)\n",
"#print(chi2s)\n",
"plt.hist(chi2s, bins=50)\n",
"print(np.mean(chi2s))"
]
},
{
"cell_type": "markdown",
"id": "ccbf4178-bb9e-440e-bed3-9365d2cee9f7",
"metadata": {},
"source": [
"Nochmal für das erste Toy"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "683060cb-9968-4bab-b512-b37621536ea2",
"metadata": {},
"outputs": [],
"source": [
"sum_xy = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*ys)))[0]\n",
"sum_y = (np.linalg.inv(J.T@W@J)@(J.T@W@(ys)))[0]\n",
"sum_x = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs)))[0]\n",
"sum_x2 = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*xs)))[0]\n",
"print(sum_xy, sum_y, sum_x, sum_x2)\n",
"mhat = (sum_xy - sum_x*sum_y)/(sum_x2-sum_x*sum_x)\n",
"ahat = (sum_y*sum_x2 - sum_x*sum_xy)/(sum_x2-sum_x*sum_x)\n",
"print(mhat, ahat)\n",
"sum_siginv = sigma_y*sigma_y/len(ys)\n",
"V_am = np.array([[sum_siginv/(sum_x2-sum_x*sum_x), -sum_siginv*sum_x/(sum_x2-sum_x*sum_x)],\n",
" [-sum_siginv*sum_x/(sum_x2-sum_x*sum_x), sum_siginv*sum_x2/(sum_x2-sum_x*sum_x)]])\n",
"print(\"chi2:\", np.sum(((ys - (xs * mhat + ahat))/sigma_y)**2))"
]
},
{
"cell_type": "markdown",
"id": "da5898ca-e7fe-4562-bed1-d9adc5866d1d",
"metadata": {},
"source": [
"# $m$ und $a$ auf Fehlerellipse \n",
"\n",
"$v = (m, a)$; $(v-\\hat v)V_{am}^{-1}(v - \\hat v)^T=1$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "73d21a5e-fdae-4e46-ab2b-c52056756f2c",
"metadata": {},
"outputs": [],
"source": [
"import scipy.optimize as opti\n",
"\n",
"sum_xy = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*ys)))[0]\n",
"sum_y = (np.linalg.inv(J.T@W@J)@(J.T@W@(ys)))[0]\n",
"sum_x = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs)))[0]\n",
"sum_x2 = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*xs)))[0]\n",
"print(sum_xy, sum_y, sum_x, sum_x2)\n",
"mhat = (sum_xy - sum_x*sum_y)/(sum_x2-sum_x*sum_x)\n",
"ahat = (sum_y*sum_x2 - sum_x*sum_xy)/(sum_x2-sum_x*sum_x)\n",
"print(mhat, ahat)\n",
"sum_siginv = sigma_y*sigma_y/len(ys)\n",
"V_am = np.array([[sum_siginv/(sum_x2-sum_x*sum_x), -sum_siginv*sum_x/(sum_x2-sum_x*sum_x)],\n",
" [-sum_siginv*sum_x/(sum_x2-sum_x*sum_x), sum_siginv*sum_x2/(sum_x2-sum_x*sum_x)]])\n",
"\n",
"def a_from_m(m, ahat, mhat, V):\n",
" Vinv = np.linalg.inv(V)\n",
" a = opti.newton(lambda a: (np.array([m - mhat, a - ahat])@Vinv@np.array([[m - mhat], [a - ahat]]))[0]-1,ahat+2)\n",
" return a\n",
"\n",
"def a_from_m2(m, ahat, mhat, V):\n",
" Vinv = np.linalg.inv(V)\n",
" a = opti.newton(lambda a: (np.array([m - mhat, a - ahat])@Vinv@np.array([[m - mhat], [a - ahat]]))[0]-1,ahat-2)\n",
" return a\n",
"\n",
"a_from_m(mhat, ahat, mhat, V_am)\n",
"m_axis = np.linspace(mhat-np.sqrt(V_am[0,0]),mhat+np.sqrt(V_am[0,0]),100)\n",
"plt.plot(m_axis, np.vectorize(a_from_m, excluded=[3])(m_axis, ahat, mhat, V_am))\n",
"plt.plot(m_axis, np.vectorize(a_from_m2, excluded=[3])(m_axis, ahat, mhat, V_am))\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "e47ddf1f-4e91-4479-8fe8-190c1eeb4628",
"metadata": {},
"source": [
"$\\chi^2$ Werte auf Ellipse"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bf1b4f5c-2634-4838-bea8-9b96e5868ffe",
"metadata": {},
"outputs": [],
"source": [
"chi2s_on_cont = np.zeros(2*len(m_axis))\n",
"for i in range(len(m_axis)):\n",
" m = m_axis[i]\n",
" a = a_from_m(m, ahat, mhat, V_am)\n",
" chi2s_on_cont[2*i] = np.sum(((ys - (xs * m + a))/sigma_y)**2)\n",
" a = a_from_m2(m, ahat, mhat, V_am)\n",
" chi2s_on_cont[2*i+1] = np.sum(((ys - (xs * m + a))/sigma_y)**2)\n",
" #print(a, m , chi2s_on_cont[2*i], chi2s_on_cont[2*i+1])\n",
"\n",
" \n",
"plt.hist(chi2s_on_cont,range=(0,10), bins=50)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "109c2f87-558b-4038-a929-4e18787a86f7",
"metadata": {},
"source": [
"$\\Delta \\chi^2$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3e43269b-a364-4a1e-9016-4fe7ef983848",
"metadata": {},
"outputs": [],
"source": [
"delta_chi2 = chi2s_on_cont - np.sum(((ys - (xs * mhat + ahat))/sigma_y)**2)\n",
"plt.hist(delta_chi2,range=(-2,2), bins=50)\n",
"print(np.mean(delta_chi2))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ce34da07-a23e-4773-b11f-4081dd3c80e0",
"metadata": {},
"outputs": [],
"source": [
"a_mesh, m_mesh = np.meshgrid(np.linspace(ahat-2*np.sqrt(V_am[1,1]),ahat+2*np.sqrt(V_am[1,1]),100), \n",
" np.linspace(mhat-2*np.sqrt(V_am[0,0]),mhat+2*np.sqrt(V_am[0,0]),100))\n",
"plt.contourf(m_mesh, a_mesh, chi2_vect(xs, ys, sigma_y, a_mesh, m_mesh)-np.sum(((ys - (xs * mhat + ahat))/sigma_y)**2),[0,1,2,3,4,5,6,7,8], cmap=cm.coolwarm)\n",
"plt.xlabel(\"m\")\n",
"plt.ylabel(\"a\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "e5a79c54-ddb4-45de-9831-f52dcae0676e",
"metadata": {},
"source": [
"# Aus scipy.optimize"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2bff85d9-ffc9-4c5e-819e-33f3ce9c56f8",
"metadata": {},
"outputs": [],
"source": [
"def fitf(x, m , a):\n",
" return m*x + a\n",
"pfit, Vfit = opti.curve_fit(fitf , xs, ys, sigma=[sigma_y]*len(ys))\n",
"print(pfit, Vfit)\n",
"#print(mhat, ahat, V_am)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "27eb4ba5-11c3-4536-8e54-450cf0814871",
"metadata": {},
"outputs": [],
"source": [
"print(Vfit/V_am)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6c235d04-988c-4a28-bb60-96d88fa3dcf7",
"metadata": {},
"outputs": [],
"source": [
"print(\"chi2 scaled:\", np.sum(((ys - (xs * pfit[0] + pfit[1]))/(sigma_y*np.sqrt(Vfit[0][0]/V_am[0][0])))**2))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c22f8191-5d24-417d-8690-c27edcb6de61",
"metadata": {},
"outputs": [],
"source": [
"pfit, Vfit = opti.curve_fit(fitf , xs, ys, sigma=[sigma_y]*len(ys), absolute_sigma=True)\n",
"print(pfit, Vfit)\n",
"#print(mhat, ahat, V_am)"
]
},
{
"cell_type": "markdown",
"id": "0e091bf2-fc2c-4bef-93af-506ad38343e5",
"metadata": {},
"source": [
"Least-squares"
]
},
{
"cell_type": "markdown",
"id": "aba1dc09-d2ce-4615-8944-9b46b69907db",
"metadata": {},
"source": [
"Minimize $Chi^2$:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "28332217-8875-49a0-be90-c7a6d00af311",
"metadata": {},
"outputs": [],
"source": [
"opti.minimize?"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6fba2b83-6cd9-4525-a15e-cb2994b1512a",
"metadata": {},
"outputs": [],
"source": [
"res = opti.minimize( lambda p: chi2(xs, ys, sigma_y, p[1], p[0]),x0=np.zeros(2))\n",
"print(res)\n",
"print(res.x, res.hess_inv)\n",
"#print(mhat, ahat, V_am)\n",
"#print(res.hess_inv/V_am)"
]
},
{
"cell_type": "markdown",
"id": "9cadcf6e-a776-4ded-bb83-23e94e03add4",
"metadata": {},
"source": [
"## Maximum-Likelihood\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d02a8b77-5f51-4a5e-bf85-2522297e36e7",
"metadata": {},
"outputs": [],
"source": [
"import scipy.stats as stats\n",
"\n",
"def like(m, a, xs, ys, sigma_y):\n",
" return np.prod(stats.norm.pdf(ys,xs*m+a,sigma_y))\n",
"\n",
"like_vect = np.vectorize(like, excluded=[2,3])\n",
"m_axis = np.linspace(-0,3,200)\n",
"\n",
"plt.plot(m_axis, like_vect(m_axis,1,xs, ys, sigma_y))\n",
"plt.xlabel(\"m\")\n",
"#plt.yscale('log')\n",
"plt.show()\n",
"plt.plot(m_axis, like_vect(2, m_axis, xs, ys, sigma_y))\n",
"plt.xlabel(\"a\")\n",
"plt.ylabel(\"$L(a; m=2)$\")\n",
"plt.savefig(\"like_a.png\")\n",
"#plt.yscale('log')\n",
"plt.show()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "09a251ef-627e-4811-85fa-dcd26ba7b589",
"metadata": {},
"outputs": [],
"source": [
"def logL(m, a, xs, ys, sigma_y):\n",
" return np.sum(-stats.norm.logpdf(ys,xs*m+a,sigma_y))\n",
"\n",
"logL_vect = np.vectorize(logL, excluded=[2,3])\n",
"\n",
"m_axis = np.linspace(-5,5,200)\n",
"plt.plot(m_axis, logL_vect(m_axis,1,xs, ys, sigma_y))\n",
"plt.xlabel(\"m\")\n",
"plt.show()\n",
"plt.plot(m_axis, logL_vect(2, m_axis, xs, ys, sigma_y))\n",
"plt.xlabel(\"a\")\n",
"plt.ylabel(\"$-\\ln L(a; m=2)$\")\n",
"plt.savefig(\"loglike_a.png\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9dc61515-4f0c-4951-a318-cb86508622e7",
"metadata": {},
"outputs": [],
"source": [
"res = opti.minimize( lambda p: logL(p[0], p[1],xs, ys, sigma_y), x0=np.zeros(2))\n",
"print(res)\n",
"print(res.x, res.hess_inv)\n",
"#print(mhat, ahat, V_am)\n",
"#print(res.hess_inv/V_am)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fa7a2440-b9f1-417e-8f5c-eb122cf9d786",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
<?xml version="1.0" encoding="utf-8" standalone="yes"?>
<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" version="1.1"
width="390px" height="390px" viewBox="-18 -120 135 135">
<title>Cumulative distribution function of log-normal distribution</title>
<defs>
<pattern id="grid" width="20" height="20" patternUnits="userSpaceOnUse">
<path d="M0,20v-20h20v20h-20" fill="none" stroke="#555" stroke-width=".25"/>
</pattern>
</defs>
<rect y="-120" width="130" height="120" fill="url(#grid)"/>
<g font-size="6" stroke="none" fill="black" text-anchor="middle" font-family="New Century Schoolbook">
<text y="6"><tspan x="0">0</tspan><tspan x="20">0.5</tspan><tspan x="40">1.0</tspan><tspan x="60">1.5</tspan><tspan x="80">2.0</tspan><tspan x="100">2.5</tspan></text>
<text x="-5" y="0">0</text><text x="-5" y="-18">0.2</text><text x="-5" y="-38">0.4</text><text x="-5" y="-58">0.6</text><text x="-5" y="-78">0.8</text><text x="-5" y="-98">1.0</text>
<g font-size="8"><text font-style="italic" x="60" y="12">x</text><text x="-12" y="-60" transform="rotate(-90,-12,-60)">CDF</text></g>
<text fill="#f00" x="105" y="-47">σ=10</text>
<text fill="#0f0" x="105" y="-67">σ=1.5</text>
<text fill="#00f" x="105" y="-77">σ=1</text>
<text fill="#aaf" x="105" y="-90">σ=0.5</text>
<text fill="#005" x="105" y="-102">σ=0.125</text>
<text fill="#070" x="60" y="-102">σ=0.25</text>
</g>
<g fill="none">
<path d="M0,0h130 M0,-120v120" stroke="black" stroke-width=".75"/>
<path d="M0,0 L0.4,-32 Q 0.47,-33 0.53,-33 0.62,-34 0.71,-34 0.82,-35 0.94,-35 1.1,-36 1.3,-36 1.5,-37 1.7,-38 1.9,-38 2.2,-39 2.6,-39 2.9,-40 3.4,-40 3.9,-41 4.6,-41 5.2,-42 6.1,-43 6.9,-43 8.1,-44 9.2,-44 11,-45 12,-45 14,-46 16,-46 19,-47 22,-48 25,-48 29,-49 34,-49 38,-50 45,-50 51,-51 59,-52 68,-52 79,-53 90,-53 110,-54 120,-54" stroke="#f00"/>
<path d="M0.4,-0.11 Q 0.45,-0.14 0.5,-0.18 0.57,-0.22 0.63,-0.28 0.71,-0.36 0.79,-0.45 0.89,-0.56 1,-0.69 1.1,-0.86 1.3,-1 1.4,-1.3 1.6,-1.5 1.8,-1.9 2,-2.2 2.2,-2.7 2.5,-3.2 2.8,-3.8 3.1,-4.4 3.5,-5.2 3.9,-6.1 4.4,-7.1 4.9,-8.1 5.6,-9.4 6.2,-11 7,-12 7.8,-14 8.8,-16 9.8,-17 11,-20 12,-22 14,-24 15,-26 17,-29 19,-31 22,-34 24,-37 27,-40 31,-43 34,-46 38,-49 43,-52 48,-55 54,-58 61,-61 68,-64 76,-67 86,-70 96,-72 110,-75 120,-77" stroke="#0f0"/>
<path d="M0,0 L 1.7,-0.074 Q1.9,-0.12 2.2,-0.19 2.6,-0.29 2.9,-0.45 3.4,-0.67 3.9,-1 4.6,-1.5 5.2,-2.1 6.1,-2.9 6.9,-4 8.1,-5.4 9.2,-7.1 11,-9.4 12,-12 14,-15 16,-18 19,-23 22,-27 25,-32 29,-37 34,-43 38,-48 45,-55 51,-60 59,-66 68,-70 79,-76 90,-79 110,-84 120,-86" stroke="#00f"/>
<path d="M0,0 H6.9 Q8.1,-0.042 9.2,-0.17 11,-0.32 12,-0.9 14,-1.7 16,-3.6 19,-6.3 22,-11 25,-17 29,-26 34,-37 38,-47 45,-60 51,-69 59,-80 68,-85 79,-93 90,-95 110,-98 120,-99" stroke="#aaf"/>
<path d="M0,0 H16 Q 19,0.077 22,-0.71 25,-1.5 29,-9.5 34,-22 38,-43 45,-71 51,-83 59,-98 68,-98 79,-100 90,-100 H120" stroke="#070"/>
<path d="M0,0 H24 Q27,0.43 30,-1.3 33,-2.7 36,-22 39,-44 42,-67 45,-87 48,-93 51,-99 54,-99 57,-100 60,-100 H120" stroke="#005"/>
</g>
</svg>
\ No newline at end of file
This diff is collapsed.
div.myheader {
position: absolute;
margin: 30px;
left: 8%;
background: blue;
font-size: xx-large;
}
div.myfooter { div.myfooter {
height: 7%;
position: absolute; position: absolute;
bottom: 0;
width: 100%;
background: red; background: red;
font-size: 120%; font-size: 150%;
right: 10%; left: 0;
} }
.rise-enabled .cm-editor { .rise-enabled .cm-editor {
......
This diff is collapsed.
%% Cell type:code id:60c7547b-61ea-4709-bb5d-30ad3c94c340 tags:
``` python
```
%% Cell type:markdown id:d2193805-8359-4380-82da-4b1cb5358290 tags: %% Cell type:markdown id:d2193805-8359-4380-82da-4b1cb5358290 tags:
# Lecture 1 # Lecture 1
--- ---
## Basic statistics ## Basic statistics
<br> <br>
<br> <br>
Hartmut Stadie Hartmut Stadie
hartmut.stadie@uni-hamburg.de hartmut.stadie@uni-hamburg.de
%% Cell type:markdown id:668ad170-00ca-4d67-ba06-71aebd7092a3 tags: %% Cell type:markdown id:668ad170-00ca-4d67-ba06-71aebd7092a3 tags:
# Parameterschätzung # Parameterschätzung
## Einführung ## Einführung
### Schätzer ### Schätzer
Schätzer $\hat a$ für $a$ aus Stichprobe $x_1,\dots x_n$ Schätzer $\hat a$ für $a$ aus Stichprobe $x_1,\dots x_n$
Anforderungen: Anforderungen:
- erwartungstreu: - erwartungstreu:
$E[\hat a]= a$ $E[\hat a]= a$
- konsistent: - konsistent:
$\lim_{n\to \infty } \hat a = a$ $\lim_{n\to \infty } \hat a = a$
- effizient: $V[\hat a]$ möglichst klein - effizient: $V[\hat a]$ möglichst klein
Schätzer für den Mittelwert: Schätzer für den Mittelwert:
$$\hat \mu = \bar x = \frac{1}{n}\sum_1^n x_i \text{ mit } V[\hat \mu] = \frac{\sigma_x^2}{N}$$ $$\hat \mu = \bar x = \frac{1}{n}\sum_1^n x_i \text{ mit } V[\hat \mu] = \frac{\sigma_x^2}{N}$$
%% Cell type:markdown id:7c7c6ecd-1611-4787-8985-da0881e90d9f tags: %% Cell type:markdown id:7c7c6ecd-1611-4787-8985-da0881e90d9f tags:
# Methode der kleinsten Quadrate # Methode der kleinsten Quadrate
## Herleitung ## Herleitung
### Methode der kleinsten Quadrate ### Methode der kleinsten Quadrate
$y(x) = mx + a$: Finde $\hat m$ und $\hat a$! $y(x) = mx + a$: Finde $\hat m$ und $\hat a$!
%% Cell type:markdown id:73e6e46a-4e41-4e66-857b-2946946498d1 tags: %% Cell type:markdown id:73e6e46a-4e41-4e66-857b-2946946498d1 tags:
<img src="./figures/11/line.png" style="width:90.0%" alt="image" /> <img src="./figures/11/line.png" style="width:90.0%" alt="image" />
%% Cell type:code id:55bf0cad-8a50-4831-8a37-daf9bcb39b71 tags: %% Cell type:code id:55bf0cad-8a50-4831-8a37-daf9bcb39b71 tags:
``` python ``` python
#hideme #hideme
import numpy as np import numpy as np
import scipy.stats as stats import scipy.stats as stats
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
def f(x): def f(x):
return 2*x + 1 return 2*x + 1
n = 10 n = 10
xs = np.linspace(0,4,n) xs = np.linspace(0,4,n)
sigma_y=0.4 sigma_y=0.4
ys = stats.multivariate_normal.rvs(f(xs), np.eye(n)*sigma_y**2, 1, random_state=42) ys = stats.multivariate_normal.rvs(f(xs), np.eye(n)*sigma_y**2, 1, random_state=42)
x_axis = np.linspace(0,4,100) x_axis = np.linspace(0,4,100)
plt.errorbar(xs,ys,yerr=sigma_y,fmt=".") plt.errorbar(xs,ys,yerr=sigma_y,fmt=".")
plt.plot(x_axis, f(x_axis),'--') plt.plot(x_axis, f(x_axis),'--')
plt.xlabel("x") plt.xlabel("x")
plt.ylabel("y") plt.ylabel("y")
plt.savefig("line.png") plt.savefig("line.png")
plt.show() plt.show()
``` ```
%% Output %% Output
%% Cell type:markdown id:64e67452-e6bd-442c-a174-e12bdb18dba0 tags: %% Cell type:markdown id:64e67452-e6bd-442c-a174-e12bdb18dba0 tags:
### Methode der kleinsten Quadrate ### Methode der kleinsten Quadrate
$$\chi^2 = \sum_i \left(\frac{y_i - \hat y(x)}{\sigma_i}\right)^2$$ $$\chi^2 = \sum_i \left(\frac{y_i - \hat y(x)}{\sigma_i}\right)^2$$
quantifiziert die Übereinstimmung von Modell zu Daten quantifiziert die Übereinstimmung von Modell zu Daten
$\rightarrow$ $\hat m$ und $\hat a$ sollten $\chi^2$ minimieren. $\rightarrow$ $\hat m$ und $\hat a$ sollten $\chi^2$ minimieren.
<img src="./figures/11/line.png" style="width:80.0%" <img src="./figures/11/line.png" style="width:80.0%"
alt="image" /> alt="image" />
%% Cell type:markdown id:bdc0949d-791e-4e29-9adc-169578e62821 tags: %% Cell type:markdown id:bdc0949d-791e-4e29-9adc-169578e62821 tags:
### Methode der kleinsten Quadrate II ### Methode der kleinsten Quadrate II
Minimiere Minimiere
$\chi^2 = \sum_i \left(\frac{y_i - \hat y(x)}{\sigma_i}\right)^2 = \sum_i \frac{(y_i - m x_i - a)^2}{\sigma_i^2}$: $\chi^2 = \sum_i \left(\frac{y_i - \hat y(x)}{\sigma_i}\right)^2 = \sum_i \frac{(y_i - m x_i - a)^2}{\sigma_i^2}$:
Erste Ableitung ist Null: Erste Ableitung ist Null:
$$\begin{aligned} $$\begin{aligned}
\frac{d\chi^2}{dm} &=& -2\sum_i x_i\frac {y_i -\hat m x_i - \hat a}{\sigma_i^2} = 0\\ \frac{d\chi^2}{dm} &=& -2\sum_i x_i\frac {y_i -\hat m x_i - \hat a}{\sigma_i^2} = 0\\
\frac{d\chi^2}{da} &=& -2\sum_i \frac{y_i - \hat m x_i - \hat a}{\sigma_i^2} = 0 \\ \frac{d\chi^2}{da} &=& -2\sum_i \frac{y_i - \hat m x_i - \hat a}{\sigma_i^2} = 0 \\
\sum_i\frac{x_iy_i}{\sigma_i^2} - \hat m \sum_i\frac{x_i^2}{\sigma_i^2}- \hat a \sum_i \frac{x_i}{\sigma_i^2} &=& 0 \\ \sum_i\frac{x_iy_i}{\sigma_i^2} - \hat m \sum_i\frac{x_i^2}{\sigma_i^2}- \hat a \sum_i \frac{x_i}{\sigma_i^2} &=& 0 \\
\sum_i\frac{y_i}{\sigma_i^2} - \hat m \sum_i\frac{x_i}{\sigma_i^2}- \hat a \sum_i \frac{1}{\sigma_i^2} &=& 0 \sum_i\frac{y_i}{\sigma_i^2} - \hat m \sum_i\frac{x_i}{\sigma_i^2}- \hat a \sum_i \frac{1}{\sigma_i^2} &=& 0
\end{aligned}$$ \end{aligned}$$
%% Cell type:markdown id:322c67fc-ddd1-4830-a5f5-6f118ef54c4c tags: %% Cell type:markdown id:322c67fc-ddd1-4830-a5f5-6f118ef54c4c tags:
### Methode der kleinsten Quadrate III ### Methode der kleinsten Quadrate III
Minimiere Minimiere
$\chi^2 = \sum_i \left(\frac{y_i - \hat y(x)}{\sigma_i}\right)^2 = \sum_i \frac{(y_i - m x_i - a)^2}{\sigma_i^2}$: $\chi^2 = \sum_i \left(\frac{y_i - \hat y(x)}{\sigma_i}\right)^2 = \sum_i \frac{(y_i - m x_i - a)^2}{\sigma_i^2}$:
$$\begin{aligned} $$\begin{aligned}
\sum_i\frac{x_iy_i}{\sigma_i^2} - \hat m \sum_i\frac{x_i^2}{\sigma_i^2}- \hat a \sum_i \frac{x_i}{\sigma_i^2} &=& 0 \\ \sum_i\frac{x_iy_i}{\sigma_i^2} - \hat m \sum_i\frac{x_i^2}{\sigma_i^2}- \hat a \sum_i \frac{x_i}{\sigma_i^2} &=& 0 \\
\sum_i\frac{y_i}{\sigma_i^2} - \hat m \sum_i\frac{x_i}{\sigma_i^2}- \hat a \sum_i \frac{1}{\sigma_i^2} &=& 0 \sum_i\frac{y_i}{\sigma_i^2} - \hat m \sum_i\frac{x_i}{\sigma_i^2}- \hat a \sum_i \frac{1}{\sigma_i^2} &=& 0
\end{aligned}$$ mit \end{aligned}$$ mit
$\frac{1}{\sum_i 1/\sigma_i^2} \sum_i \frac{f}{\sigma_i^2} = \langle f \rangle$: $\frac{1}{\sum_i 1/\sigma_i^2} \sum_i \frac{f}{\sigma_i^2} = \langle f \rangle$:
$$\begin{aligned} $$\begin{aligned}
\langle xy \rangle -\langle x^2 \rangle \hat m& - \langle x \rangle \hat a&= 0\\ \langle xy \rangle -\langle x^2 \rangle \hat m& - \langle x \rangle \hat a&= 0\\
\langle y \rangle - \langle x \rangle \hat m& - \hat a& = 0 \langle y \rangle - \langle x \rangle \hat m& - \hat a& = 0
\end{aligned}$$ \end{aligned}$$
%% Cell type:markdown id:fb25687c-4410-4281-b540-39369732fb26 tags: %% Cell type:markdown id:fb25687c-4410-4281-b540-39369732fb26 tags:
### Methode der kleinsten Quadrate IV ### Methode der kleinsten Quadrate IV
$$\begin{aligned} $$\begin{aligned}
\hat m&=&\frac{\langle xy \rangle - \langle y \rangle\langle x \rangle}{\langle x^2 \rangle - \langle x \rangle^2} = \frac{1}{\sum_i 1/\sigma_i^2} \sum_i \frac{x_i - \langle x \rangle}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)}y_i\\ \hat m&=&\frac{\langle xy \rangle - \langle y \rangle\langle x \rangle}{\langle x^2 \rangle - \langle x \rangle^2} = \frac{1}{\sum_i 1/\sigma_i^2} \sum_i \frac{x_i - \langle x \rangle}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)}y_i\\
\hat a &=& \frac{ \langle y \rangle \langle x^2 \rangle- \langle y \rangle \langle x \rangle^2- \langle x \rangle \langle xy \rangle+ \langle y \rangle \langle x \rangle^2}{ \langle x^2 \rangle- \langle x \rangle^2}\\ \hat a &=& \frac{ \langle y \rangle \langle x^2 \rangle- \langle y \rangle \langle x \rangle^2- \langle x \rangle \langle xy \rangle+ \langle y \rangle \langle x \rangle^2}{ \langle x^2 \rangle- \langle x \rangle^2}\\
&=& \frac{ \langle y \rangle \langle x^2 \rangle - \langle x \rangle \langle xy \rangle}{ \langle x^2 \rangle - \langle x \rangle^2} = \frac{1}{\sum_i 1/\sigma_i^2} \sum_i \frac{\langle x^2 \rangle - \langle x \rangle x_i}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)}y_i &=& \frac{ \langle y \rangle \langle x^2 \rangle - \langle x \rangle \langle xy \rangle}{ \langle x^2 \rangle - \langle x \rangle^2} = \frac{1}{\sum_i 1/\sigma_i^2} \sum_i \frac{\langle x^2 \rangle - \langle x \rangle x_i}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)}y_i
\end{aligned}$$ \end{aligned}$$
%% Cell type:markdown id:b08ade05-e5eb-4a7e-9315-31a45c51aadb tags: %% Cell type:markdown id:b08ade05-e5eb-4a7e-9315-31a45c51aadb tags:
## Fehler ## Fehler
### Fehler ### Fehler
$$\begin{aligned} $$\begin{aligned}
V(\hat m) = \sum_i \left(\frac{d\hat m}{y_i}\sigma_i\right)^2\text{; }\frac{d\hat m}{y_i} & = & \frac{1}{\sum_i 1/\sigma_i^2} \frac{x_i - \langle x \rangle}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)} \\ V(\hat m) = \sum_i \left(\frac{d\hat m}{y_i}\sigma_i\right)^2\text{; }\frac{d\hat m}{y_i} & = & \frac{1}{\sum_i 1/\sigma_i^2} \frac{x_i - \langle x \rangle}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)} \\
V(\hat a) = \sum_i \left(\frac{d\hat a}{y_i}\sigma_i\right)^2\text{; }\frac{d\hat a}{y_i} & = & \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle - \langle x \rangle x_i}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)} V(\hat a) = \sum_i \left(\frac{d\hat a}{y_i}\sigma_i\right)^2\text{; }\frac{d\hat a}{y_i} & = & \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle - \langle x \rangle x_i}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)}
\end{aligned}$$ $$\begin{aligned} \end{aligned}$$ $$\begin{aligned}
V(\hat m) &=& \left(\frac{1}{\sum_i 1/\sigma_i^2}\right)^2 \sum_i \left(\frac{x_i - \langle x \rangle}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)}\right)^2 \sigma_i^2 \\ V(\hat m) &=& \left(\frac{1}{\sum_i 1/\sigma_i^2}\right)^2 \sum_i \left(\frac{x_i - \langle x \rangle}{\sigma_i^2(\langle x^2 \rangle - \langle x \rangle^2)}\right)^2 \sigma_i^2 \\
&=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle - 2\langle x \rangle \langle x \rangle + \langle x \rangle^2}{(\langle x^2 \rangle - \langle x \rangle^2)^2} &=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle - 2\langle x \rangle \langle x \rangle + \langle x \rangle^2}{(\langle x^2 \rangle - \langle x \rangle^2)^2}
= \frac{1}{\sum_i 1/\sigma_i^2} \frac{1}{\langle x^2 \rangle - \langle x \rangle^2} \\ = \frac{1}{\sum_i 1/\sigma_i^2} \frac{1}{\langle x^2 \rangle - \langle x \rangle^2} \\
V(\hat a) &=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle^2 - 2\langle x^2 \rangle\langle x \rangle^2 + \langle x^2 \rangle\langle x \rangle^2}{(\langle x^2 \rangle - \langle x \rangle^2)^2} V(\hat a) &=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle^2 - 2\langle x^2 \rangle\langle x \rangle^2 + \langle x^2 \rangle\langle x \rangle^2}{(\langle x^2 \rangle - \langle x \rangle^2)^2}
= \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle}{\langle x^2 \rangle - \langle x \rangle^2} = \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle}{\langle x^2 \rangle - \langle x \rangle^2}
\end{aligned}$$ \end{aligned}$$
%% Cell type:markdown id:d8828d04-0af8-4dc3-a846-24acbc9a0f8e tags: %% Cell type:markdown id:d8828d04-0af8-4dc3-a846-24acbc9a0f8e tags:
### Korrelation ### Korrelation
$$\begin{aligned} $$\begin{aligned}
V(\hat m) &=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{1}{\langle x^2 \rangle - \langle x \rangle^2} \\ V(\hat m) &=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{1}{\langle x^2 \rangle - \langle x \rangle^2} \\
V(\hat a) &=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle}{\langle x^2 \rangle - \langle x \rangle^2}\\ V(\hat a) &=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle}{\langle x^2 \rangle - \langle x \rangle^2}\\
\text{cov}(\hat m, \hat a) &=&= \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle (x-\langle x \rangle)(\langle x^2 \rangle - \langle x \rangle x)\rangle}{(\langle x^2 \rangle - \langle x \rangle^2)^2}\\ \text{cov}(\hat m, \hat a) &=&= \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle (x-\langle x \rangle)(\langle x^2 \rangle - \langle x \rangle x)\rangle}{(\langle x^2 \rangle - \langle x \rangle^2)^2}\\
&=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle \langle x \rangle - \langle x \rangle \langle x^2 \rangle - \langle x \rangle \langle x^2 \rangle + \langle x \rangle^2\langle x \rangle}{(\langle x^2 \rangle - \langle x \rangle^2)^2}\\ &=& \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x^2 \rangle \langle x \rangle - \langle x \rangle \langle x^2 \rangle - \langle x \rangle \langle x^2 \rangle + \langle x \rangle^2\langle x \rangle}{(\langle x^2 \rangle - \langle x \rangle^2)^2}\\
&=& - \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x \rangle}{\langle x^2 \rangle - \langle x \rangle^2} &=& - \frac{1}{\sum_i 1/\sigma_i^2} \frac{\langle x \rangle}{\langle x^2 \rangle - \langle x \rangle^2}
\end{aligned}$$ \end{aligned}$$
### Beispiel in Jupyter ### Beispiel in Jupyter
%% Cell type:markdown id:424bdd1f-53bf-422b-bc4a-0702231b976d tags: %% Cell type:markdown id:424bdd1f-53bf-422b-bc4a-0702231b976d tags:
### Minimales $\chi^2$ ### Minimales $\chi^2$
$$\begin{aligned} $$\begin{aligned}
\chi^2 &=& \sum_i \frac{(y_i - \hat m x_i - \hat a)^2}{\sigma_i^2} = \sum_i \frac{\left[y_i - \frac{\langle xy \rangle - \langle y \rangle\langle x \rangle}{\langle x^2 \rangle - \langle x \rangle^2} x_i - \frac{ \langle y \rangle \langle x^2 \rangle - \langle x \rangle \langle xy \rangle}{ \langle x^2 \rangle - \langle x \rangle^2} \right]^2}{\sigma_i^2}\\ \chi^2 &=& \sum_i \frac{(y_i - \hat m x_i - \hat a)^2}{\sigma_i^2} = \sum_i \frac{\left[y_i - \frac{\langle xy \rangle - \langle y \rangle\langle x \rangle}{\langle x^2 \rangle - \langle x \rangle^2} x_i - \frac{ \langle y \rangle \langle x^2 \rangle - \langle x \rangle \langle xy \rangle}{ \langle x^2 \rangle - \langle x \rangle^2} \right]^2}{\sigma_i^2}\\
& = & \sum_i \frac{\left[(\langle x^2 \rangle - \langle x \rangle^2)y_i - (\langle xy \rangle - \langle y \rangle\langle x \rangle)x_i - \langle y \rangle \langle x^2 \rangle + \langle x \rangle \langle xy \rangle\right]^2}{\sigma_i^2 ( \langle x^2 \rangle - \langle x \rangle^2)^2} \\ & = & \sum_i \frac{\left[(\langle x^2 \rangle - \langle x \rangle^2)y_i - (\langle xy \rangle - \langle y \rangle\langle x \rangle)x_i - \langle y \rangle \langle x^2 \rangle + \langle x \rangle \langle xy \rangle\right]^2}{\sigma_i^2 ( \langle x^2 \rangle - \langle x \rangle^2)^2} \\
&=& \dots\\ &=& \dots\\
& =& (\sum_i \frac{1}{\sigma_i^2}) V(y) ( 1- \rho^2_{xy}) & =& (\sum_i \frac{1}{\sigma_i^2}) V(y) ( 1- \rho^2_{xy})
\end{aligned}$$ \end{aligned}$$
%% Cell type:markdown id:e89f415e-8836-4b97-893c-a7335c3a21e5 tags: %% Cell type:markdown id:e89f415e-8836-4b97-893c-a7335c3a21e5 tags:
### Beispiel in Jupyter ### Beispiel in Jupyter
## In Python ## In Python
%% Cell type:markdown id:7311c0ff-0ce0-4427-a50e-b6698100e454 tags: %% Cell type:markdown id:7311c0ff-0ce0-4427-a50e-b6698100e454 tags:
### Mit Python I ### Mit Python I
Mit scipy.optimize: Mit scipy.optimize:
``` ```
import scipy.optimize as opti import scipy.optimize as opti
def fitf(x, m , a): def fitf(x, m , a):
return m*x + a return m*x + a
pfit, Vfit = opti.curve_fit(fitf , xs, ys, pfit, Vfit = opti.curve_fit(fitf , xs, ys,
sigma=[sigma_y]*len(ys),absolue_sigma=True) sigma=[sigma_y]*len(ys),absolue_sigma=True)
print(pfit, Vfit) print(pfit, Vfit)
``` ```
Vorsicht! Falsche Unsicherheit ohne `absolute_sigma=True` Vorsicht! Falsche Unsicherheit ohne `absolute_sigma=True`
%% Cell type:markdown id:a5fec52e-bd3a-4437-bb43-620de44939b2 tags: %% Cell type:markdown id:a5fec52e-bd3a-4437-bb43-620de44939b2 tags:
### Mit Python II ### Mit Python II
Mit scipy.optimize: Mit scipy.optimize:
``` ```
def chi2(x, y, sy, a, m): def chi2(x, y, sy, a, m):
my = m * x + a my = m * x + a
r = (y - my)/sy r = (y - my)/sy
return np.sum(r**2) return np.sum(r**2)
res = opti.minimize( lambda p: chi2(xs, ys, sigma_y, p[1], p[0]),x0=np.zeros(2)) res = opti.minimize( lambda p: chi2(xs, ys, sigma_y, p[1], p[0]),x0=np.zeros(2))
print(res.x, res.hess_inv*2) print(res.x, res.hess_inv*2)
``` ```
%% Cell type:markdown id:84b00301-c0d3-4595-858e-4f784d69c876 tags: %% Cell type:markdown id:84b00301-c0d3-4595-858e-4f784d69c876 tags:
### Inverse Hesse-Matrix und $\chi^2$ ### Inverse Hesse-Matrix und $\chi^2$
$\Delta \chi2$ und Kovarianz Ellipse um Minimum gemäß Kovarianzmatrix $\Delta \chi2$ und Kovarianz Ellipse um Minimum gemäß Kovarianzmatrix
genau bei $\Delta \chi^2 = 1$. genau bei $\Delta \chi^2 = 1$.
$$1 = \delta \chi^2 = (\vec a -\hat \vec a)^T V^{-1} (\vec a-\hat \vec a)$$ $$1 = \delta \chi^2 = (\vec a -\hat \vec a)^T V^{-1} (\vec a-\hat \vec a)$$
Mit Mit
$\chi^2(\vec a) = \chi^2(\hat \vec a) + (\vec a -\hat \vec a)^T V^{-1} (\vec a-\hat \vec a)$ $\chi^2(\vec a) = \chi^2(\hat \vec a) + (\vec a -\hat \vec a)^T V^{-1} (\vec a-\hat \vec a)$
und und
$H_{ij} = \frac{\partial^2 \chi^2(\vec a)}{\partial a_i \partial a_j}$ $H_{ij} = \frac{\partial^2 \chi^2(\vec a)}{\partial a_i \partial a_j}$
$$H_{ij} = \frac{\partial^2 (a_k -\hat a_k) V^{-1}_{kl} (a_l -\hat a_l)}{\partial a_i \partial a_j} = \frac{\partial( \delta_{ik}V^{-1}_{kl} (a_l -\hat a_l) + (a_k -\hat a_k) V^{-1}_{kl} \delta_{il})}{\partial a_j}$$ $$H_{ij} = \frac{\partial^2 (a_k -\hat a_k) V^{-1}_{kl} (a_l -\hat a_l)}{\partial a_i \partial a_j} = \frac{\partial( \delta_{ik}V^{-1}_{kl} (a_l -\hat a_l) + (a_k -\hat a_k) V^{-1}_{kl} \delta_{il})}{\partial a_j}$$
$$H_{ij} = \delta_{ik}V^{-1}_{kl}\delta_{lj} + \delta_{jk}V^{-1}_{kl}\delta_{il} = 2V^{-1}_{ij} \text{ und } V_{ij} = 2 * H^{-1}_{ij}$$ $$H_{ij} = \delta_{ik}V^{-1}_{kl}\delta_{lj} + \delta_{jk}V^{-1}_{kl}\delta_{il} = 2V^{-1}_{ij} \text{ und } V_{ij} = 2 * H^{-1}_{ij}$$
Vorsicht! Manche Algorithmen in `minimize` berechnen keine inverse Vorsicht! Manche Algorithmen in `minimize` berechnen keine inverse
Hesse-Matrix. Hesse-Matrix.
%% Cell type:markdown id:632583d4-1e97-498b-b8f3-91e259baa24d tags: %% Cell type:markdown id:632583d4-1e97-498b-b8f3-91e259baa24d tags:
# Maximum-Likelihood # Maximum-Likelihood
Maximum-Likelihood (ML) Daten: $x_1,...,x_N$ Maximum-Likelihood (ML) Daten: $x_1,...,x_N$
Wahrscheinlichkeit der Daten für Modell mit Parametern $a$: Wahrscheinlichkeit der Daten für Modell mit Parametern $a$:
$$P(x_1,...,x_N; a) = \prod_i P(x_i ; a)$$ $$P(x_1,...,x_N; a) = \prod_i P(x_i ; a)$$
Likelihoodfunktion: $$L(a) = \prod_i P(x_i ; a)$$ Likelihoodfunktion: $$L(a) = \prod_i P(x_i ; a)$$
ML-Schätzer $\hat a$: Maximum von $L(a)$: ML-Schätzer $\hat a$: Maximum von $L(a)$:
$$\left.\frac{dL}{da}\right|_{a = \hat a} = 0$$ (praktischer: $$\left.\frac{dL}{da}\right|_{a = \hat a} = 0$$ (praktischer:
Log-Likelihood: $-\ln L = \sum_i -\ln P(x_i; a)$) Log-Likelihood: $-\ln L = \sum_i -\ln P(x_i; a)$)
%% Cell type:markdown id:dd0afbf7-8504-46f8-b3da-4fb33487a635 tags: %% Cell type:markdown id:dd0afbf7-8504-46f8-b3da-4fb33487a635 tags:
### Beispiel ### Beispiel
$y(x) = mx + a$: Finde $\hat m$ und $\hat a$ Daten: $y_1,...,y_N$ und $y(x) = mx + a$: Finde $\hat m$ und $\hat a$ Daten: $y_1,...,y_N$ und
Modell: $$P(y_i; m, a) = G(y_i; \mu = m x_i + a, \sigma=\sigma_i)$$ Modell: $$P(y_i; m, a) = G(y_i; \mu = m x_i + a, \sigma=\sigma_i)$$
$$L(m, a) = \prod_i G(y_i; \mu = m x_i + a, \sigma=\sigma_i)$$ $$L(m, a) = \prod_i G(y_i; \mu = m x_i + a, \sigma=\sigma_i)$$
0.5 <img src="./figures/11/line.png" alt="image" /> 0.5 <img src="./figures/11/line.png" alt="image" />
<img src="./figures/11/like_a.png" style="width:49.0%" <img src="./figures/11/like_a.png" style="width:49.0%"
alt="image" /> alt="image" />
<img src="./figures/11/loglike_a.png" style="width:49.0%" <img src="./figures/11/loglike_a.png" style="width:49.0%"
alt="image" /> alt="image" />
ML-Schätzer für Poisson $\mu$ $$\begin{aligned} ML-Schätzer für Poisson $\mu$ $$\begin{aligned}
L(\mu) & = & \prod_i^N P(k_i; \mu) = \prod_i^N \frac{\mu^{k_i}e^{-\mu}}{k_i!}\\ L(\mu) & = & \prod_i^N P(k_i; \mu) = \prod_i^N \frac{\mu^{k_i}e^{-\mu}}{k_i!}\\
\ln L(\mu) & = & \sum_{i=1}^N \left( \ln \mu^{k_i} + \ln e^{-\mu} - \ln k_i!\right)\\ \ln L(\mu) & = & \sum_{i=1}^N \left( \ln \mu^{k_i} + \ln e^{-\mu} - \ln k_i!\right)\\
& = & \sum_{i=1}^N \left( k_i \ln \mu -\mu - \ln k_i!\right)\\ & = & \sum_{i=1}^N \left( k_i \ln \mu -\mu - \ln k_i!\right)\\
0 \stackrel{!}{=} \frac{d \ln L(\mu)}{d\mu} \Big|_{\hat \mu}& = & \sum_{i=1}^N \left( \frac{k_i}{\hat \mu} - 1\right) = \sum_{i=1}^N \frac{k_i}{\hat\mu} - N\\ 0 \stackrel{!}{=} \frac{d \ln L(\mu)}{d\mu} \Big|_{\hat \mu}& = & \sum_{i=1}^N \left( \frac{k_i}{\hat \mu} - 1\right) = \sum_{i=1}^N \frac{k_i}{\hat\mu} - N\\
N & = & \frac{1}{\hat\mu} \sum_{i=1}^N k_i \rightarrow \hat\mu = \frac{1} {N} \sum_{i=1}^N k_i N & = & \frac{1}{\hat\mu} \sum_{i=1}^N k_i \rightarrow \hat\mu = \frac{1} {N} \sum_{i=1}^N k_i
\end{aligned}$$ \end{aligned}$$
%% Cell type:markdown id:d9571970-772e-4e95-8474-82e0afb1dd68 tags: %% Cell type:markdown id:d9571970-772e-4e95-8474-82e0afb1dd68 tags:
Varianz des ML-Schätzers Varianz des ML-Schätzers
Rao-Cramér-Frechet-Ungleichung: Schätzer $\hat a$ mit Bias (Verzerrung) Rao-Cramér-Frechet-Ungleichung: Schätzer $\hat a$ mit Bias (Verzerrung)
$b$ $b$
$$V(\hat a) \geq \frac{\left(1+ \frac{\partial b}{\partial a} \right)^2}{E\left[-\frac{\partial^2 \ln L}{\partial a^2}\right]}$$ $$V(\hat a) \geq \frac{\left(1+ \frac{\partial b}{\partial a} \right)^2}{E\left[-\frac{\partial^2 \ln L}{\partial a^2}\right]}$$
Fisher-Information: Fisher-Information:
$$I(\hat a) = E\left [-\frac{\partial^2 \ln L}{\partial a^2}\right]$$ $$I(\hat a) = E\left [-\frac{\partial^2 \ln L}{\partial a^2}\right]$$
ML-Schätzer für Poisson $V(\hat \mu)$ $$\begin{aligned} ML-Schätzer für Poisson $V(\hat \mu)$ $$\begin{aligned}
V(\hat \mu) & \geq &\frac{\left(1+ \frac{\partial b}{\partial \mu} \right)^2}{E\left[-\frac{\partial^2 \ln L}{\mu^2}\right]} \\ V(\hat \mu) & \geq &\frac{\left(1+ \frac{\partial b}{\partial \mu} \right)^2}{E\left[-\frac{\partial^2 \ln L}{\mu^2}\right]} \\
& = & \frac{1}{E\left[-\frac{\partial(\sum_{i=1}^N \frac{k_i}{\mu} - N)}{\partial \mu^2}\right]} \\ & = & \frac{1}{E\left[-\frac{\partial(\sum_{i=1}^N \frac{k_i}{\mu} - N)}{\partial \mu^2}\right]} \\
& = & \frac{1}{E\left[-\sum_{i=1}^N \frac{-k_i}{\hat \mu^2}\right]} = \frac{1}{E\left[\sum_{i=1}^N \frac{k_i}{\hat \mu^2}\right]} \\ & = & \frac{1}{E\left[-\sum_{i=1}^N \frac{-k_i}{\hat \mu^2}\right]} = \frac{1}{E\left[\sum_{i=1}^N \frac{k_i}{\hat \mu^2}\right]} \\
& = & \frac{1}{\frac{1}{\hat \mu^2}E\left[\sum_{i=1}^N k_i \right]} = \frac{1}{\frac{1}{\hat \mu^2}E\left[N \hat \mu \right]}\\ & = & \frac{1}{\frac{1}{\hat \mu^2}E\left[\sum_{i=1}^N k_i \right]} = \frac{1}{\frac{1}{\hat \mu^2}E\left[N \hat \mu \right]}\\
& = & \frac{\hat \mu}{N} & = & \frac{\hat \mu}{N}
\end{aligned}$$ \end{aligned}$$
Varianz für mehrere Parameter $\vec \theta$ Varianz für mehrere Parameter $\vec \theta$
Für effizienten und erwartungstreuen Schätzer: Für effizienten und erwartungstreuen Schätzer:
$$\left(V^{-1}\right)_{ij} = E\left[ -\frac{\partial^2 \ln L(\theta)}{\partial \theta_i \partial \theta_j}\right]$$ $$\left(V^{-1}\right)_{ij} = E\left[ -\frac{\partial^2 \ln L(\theta)}{\partial \theta_i \partial \theta_j}\right]$$
Näherung für große Datensätze: Näherung für große Datensätze:
$$\left(\hat V^{-1}\right)_{ij} = -\frac{\partial^2 \ln L(\theta)}{\partial \theta_i \partial \theta_j}\Big|_{\theta=\hat \theta} =$$ $$\left(\hat V^{-1}\right)_{ij} = -\frac{\partial^2 \ln L(\theta)}{\partial \theta_i \partial \theta_j}\Big|_{\theta=\hat \theta} =$$
Graphisch: Graphisch:
$$\ln L(\theta) \approx \ln L(\hat \theta) + \frac{\partial \ln L}{\partial \theta}\Big|_{\hat \theta}(\theta - \hat \theta) + \frac{1}{2} \frac{\partial^2 \ln L}{\partial \theta^2}(\theta - \hat \theta)^2$$ $$\ln L(\theta) \approx \ln L(\hat \theta) + \frac{\partial \ln L}{\partial \theta}\Big|_{\hat \theta}(\theta - \hat \theta) + \frac{1}{2} \frac{\partial^2 \ln L}{\partial \theta^2}(\theta - \hat \theta)^2$$
$$\ln L(\hat \theta + \sigma_\theta) \approx \ln L(\hat \theta) + \frac{1}{2} \frac{\partial^2 \ln L}{\partial \theta^2}(\sigma_\theta)^2 = \ln L(\hat \theta) - \frac{1}{2}$$ $$\ln L(\hat \theta + \sigma_\theta) \approx \ln L(\hat \theta) + \frac{1}{2} \frac{\partial^2 \ln L}{\partial \theta^2}(\sigma_\theta)^2 = \ln L(\hat \theta) - \frac{1}{2}$$
Zusammenhang ML und $\chi^2$ Zusammenhang ML und $\chi^2$
Likelihood-Quotient: Likelihood-Quotient:
$$\lambda = -2 \ln \frac{L(\hat \theta)}{L(\hat \theta^\prime_\text{saturiert})}$$ $$\lambda = -2 \ln \frac{L(\hat \theta)}{L(\hat \theta^\prime_\text{saturiert})}$$
Mit Normalverteilung: $$\begin{aligned} Mit Normalverteilung: $$\begin{aligned}
\lambda &=& -2 \ln \frac{L(\hat \theta)}{L(\hat \theta^\prime_\text{saturiert})} = -2 \ln \frac{\prod_i G(x_i; \hat \mu, \sigma_i)}{\prod_i G(x_i; x_i, \sigma_i)}\\ \lambda &=& -2 \ln \frac{L(\hat \theta)}{L(\hat \theta^\prime_\text{saturiert})} = -2 \ln \frac{\prod_i G(x_i; \hat \mu, \sigma_i)}{\prod_i G(x_i; x_i, \sigma_i)}\\
& = & -2 \ln \frac{\frac{1}{\sqrt{2\pi}\sigma_i}exp\left(\frac{(x_i-\hat \mu)^2}{2\sigma_i^2}\right)}{\frac{1}{\sqrt{2\pi}\sigma_i}exp\left(\frac{(x_i-x_i)^2}{2\sigma_i^2}\right)} = -2\ln exp\left(\frac{(x_i-\hat \mu)^2}{2\sigma_i^2}\right) \\ & = & -2 \ln \frac{\frac{1}{\sqrt{2\pi}\sigma_i}exp\left(\frac{(x_i-\hat \mu)^2}{2\sigma_i^2}\right)}{\frac{1}{\sqrt{2\pi}\sigma_i}exp\left(\frac{(x_i-x_i)^2}{2\sigma_i^2}\right)} = -2\ln exp\left(\frac{(x_i-\hat \mu)^2}{2\sigma_i^2}\right) \\
& = & -2 \frac{(x_i-\hat \mu)^2}{2\sigma_i^2} = \chi^2 \text{; also } \ln L(\theta) = - \chi^2(\theta) / 2 & = & -2 \frac{(x_i-\hat \mu)^2}{2\sigma_i^2} = \chi^2 \text{; also } \ln L(\theta) = - \chi^2(\theta) / 2
\end{aligned}$$ \end{aligned}$$
%% Cell type:markdown id:7ae85d99-d6fa-409b-abb2-144cc24570db tags: %% Cell type:markdown id:7ae85d99-d6fa-409b-abb2-144cc24570db tags:
# Zusammenfassung und Ausblick # Zusammenfassung und Ausblick
## Zusammenfassung und Ausblick ## Zusammenfassung und Ausblick
Zusammenfassung Zusammenfassung
- Methode der kleinsten Quadrate ($\chi^2$) - Methode der kleinsten Quadrate ($\chi^2$)
- Maximum-Likelihood - Maximum-Likelihood
- Zusammenhang $\chi^2$-ML - Zusammenhang $\chi^2$-ML
- Minimierung - Minimierung
- Literatur: - Literatur:
- Glen Cowan, Statistical Data Analysis, - Glen Cowan, Statistical Data Analysis,
[pdf](https://www.sherrytowers.com/cowan_statistical_data_analysis.pdf) [pdf](https://www.sherrytowers.com/cowan_statistical_data_analysis.pdf)
- Roger John Barlow, Statistics: A Guide to the Use of Statistical - Roger John Barlow, Statistics: A Guide to the Use of Statistical
Methods in the Physical Sciences, Methods in the Physical Sciences,
[Skript](https://arxiv.org/pdf/1905.12362.pdf) [Skript](https://arxiv.org/pdf/1905.12362.pdf)
- Volker Blobel, Erich Lohrmann, Statistische und numerische - Volker Blobel, Erich Lohrmann, Statistische und numerische
Methoden der Datenanalyse, Methoden der Datenanalyse,
[pdf](https://www.desy.de/~sschmitt/blobel/eBuch.pdf) [pdf](https://www.desy.de/~sschmitt/blobel/eBuch.pdf)
# Bibliography # Bibliography
Bibliography Bibliography
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment