From a58781f3e1d3e1fbd7db993cd0a1a9862edc146c Mon Sep 17 00:00:00 2001 From: Hartmut Stadie <hartmut.stadie@cern.ch> Date: Wed, 2 Oct 2024 14:10:55 +0200 Subject: [PATCH] new requirements --- environment.yml | 847 ++++++++++++------------------------------------ 1 file changed, 213 insertions(+), 634 deletions(-) diff --git a/environment.yml b/environment.yml index 72e3694..0f7e35c 100644 --- a/environment.yml +++ b/environment.yml @@ -1,634 +1,213 @@ -{ - "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 -} +name: jupy +channels: + - conda-forge + - defaults +dependencies: + - anyio=4.4.0=pyhd8ed1ab_0 + - appnope=0.1.4=pyhd8ed1ab_0 + - argon2-cffi=23.1.0=pyhd8ed1ab_0 + - argon2-cffi-bindings=21.2.0=py39h06df861_5 + - arrow=1.3.0=pyhd8ed1ab_0 + - asttokens=2.4.1=pyhd8ed1ab_0 + - async-lru=2.0.4=pyhd8ed1ab_0 + - attrs=23.2.0=pyh71513ae_0 + - babel=2.14.0=pyhd8ed1ab_0 + - beautifulsoup4=4.12.3=pyha770c72_0 + - bleach=6.1.0=pyhd8ed1ab_0 + - brotli=1.1.0=hb547adb_1 + - brotli-bin=1.1.0=hb547adb_1 + - brotli-python=1.1.0=py39hb198ff7_1 + - bzip2=1.0.8=h93a5062_5 + - ca-certificates=2024.8.30=hf0a4a13_0 + - cached-property=1.5.2=hd8ed1ab_1 + - cached_property=1.5.2=pyha770c72_1 + - certifi=2024.8.30=pyhd8ed1ab_0 + - cffi=1.17.1=py39h7f933ea_0 + - charset-normalizer=3.3.2=pyhd8ed1ab_0 + - comm=0.2.2=pyhd8ed1ab_0 + - contourpy=1.3.0=py39h157d57c_2 + - cycler=0.12.1=pyhd8ed1ab_0 + - debugpy=1.8.6=py39hfa9831e_0 + - decorator=5.1.1=pyhd8ed1ab_0 + - defusedxml=0.7.1=pyhd8ed1ab_0 + - entrypoints=0.4=pyhd8ed1ab_0 + - exceptiongroup=1.2.0=pyhd8ed1ab_2 + - executing=2.0.1=pyhd8ed1ab_0 + - fonttools=4.54.1=py39h06df861_0 + - fqdn=1.5.1=pyhd8ed1ab_0 + - freetype=2.12.1=hadb7bae_2 + - h11=0.14.0=pyhd8ed1ab_0 + - h2=4.1.0=pyhd8ed1ab_0 + - hpack=4.0.0=pyh9f0ad1d_0 + - httpcore=1.0.5=pyhd8ed1ab_0 + - httpx=0.27.0=pyhd8ed1ab_0 + - hyperframe=6.0.1=pyhd8ed1ab_0 + - icu=75.1=hfee45f7_0 + - idna=3.7=pyhd8ed1ab_0 + - importlib-metadata=8.0.0=pyha770c72_0 + - importlib-resources=6.4.0=pyhd8ed1ab_0 + - importlib_metadata=8.0.0=hd8ed1ab_0 + - importlib_resources=6.4.0=pyhd8ed1ab_0 + - ipykernel=6.29.5=pyh57ce528_0 + - ipython=8.18.1=pyh707e725_3 + - ipython_genutils=0.2.0=pyhd8ed1ab_1 + - ipywidgets=8.1.3=pyhd8ed1ab_0 + - isoduration=20.11.0=pyhd8ed1ab_0 + - jedi=0.19.1=pyhd8ed1ab_0 + - jinja2=3.1.4=pyhd8ed1ab_0 + - json5=0.9.25=pyhd8ed1ab_0 + - jsonpointer=3.0.0=py39h2804cbe_1 + - jsonschema=4.22.0=pyhd8ed1ab_0 + - jsonschema-specifications=2023.12.1=pyhd8ed1ab_0 + - jsonschema-with-format-nongpl=4.22.0=pyhd8ed1ab_0 + - jupyter=1.1.1=pyhd8ed1ab_0 + - jupyter-lsp=2.2.5=pyhd8ed1ab_0 + - jupyter_client=7.4.9=pyhd8ed1ab_0 + - jupyter_console=6.6.3=pyhd8ed1ab_0 + - jupyter_contrib_core=0.4.0=pyhd8ed1ab_0 + - jupyter_contrib_nbextensions=0.7.0=pyhd8ed1ab_0 + - jupyter_core=5.7.2=pyh31011fe_1 + - jupyter_events=0.10.0=pyhd8ed1ab_0 + - jupyter_highlight_selected_word=0.2.0=pyhd8ed1ab_1006 + - jupyter_latex_envs=1.4.6=pyhd8ed1ab_1002 + - jupyter_nbextensions_configurator=0.6.1=pyhd8ed1ab_0 + - jupyter_server=2.14.1=pyhd8ed1ab_0 + - jupyter_server_terminals=0.5.3=pyhd8ed1ab_0 + - jupyterlab=4.2.5=pyhd8ed1ab_0 + - jupyterlab_pygments=0.3.0=pyhd8ed1ab_1 + - jupyterlab_rise=0.42.0=pyhd8ed1ab_0 + - jupyterlab_server=2.27.2=pyhd8ed1ab_0 + - jupyterlab_widgets=3.0.11=pyhd8ed1ab_0 + - kiwisolver=1.4.7=py39h157d57c_0 + - krb5=1.21.3=h237132a_0 + - lcms2=2.16=ha0e7c42_0 + - lerc=4.0.0=h9a09cb3_0 + - libblas=3.9.0=22_osxarm64_openblas + - libbrotlicommon=1.1.0=hb547adb_1 + - libbrotlidec=1.1.0=hb547adb_1 + - libbrotlienc=1.1.0=hb547adb_1 + - libcblas=3.9.0=22_osxarm64_openblas + - libcxx=17.0.6=h0812c0d_3 + - libdeflate=1.20=h93a5062_0 + - libedit=3.1.20191231=hc8eb9b7_2 + - libexpat=2.6.3=hf9b8971_0 + - libffi=3.4.2=h3422bc3_5 + - libgfortran=5.0.0=13_2_0_hd922786_3 + - libgfortran5=13.2.0=hf226fd6_3 + - libiconv=1.17=h0d3ecfb_2 + - libjpeg-turbo=3.0.0=hb547adb_1 + - liblapack=3.9.0=22_osxarm64_openblas + - libopenblas=0.3.27=openmp_h517c56d_1 + - libpng=1.6.43=h091b4b1_0 + - libsodium=1.0.18=h27ca646_1 + - libsqlite=3.46.1=hc14010f_0 + - libtiff=4.6.0=h07db509_3 + - libuv=1.48.0=h93a5062_0 + - libwebp-base=1.4.0=h93a5062_0 + - libxcb=1.16=hf2054a2_0 + - libxml2=2.12.7=h01dff8b_4 + - libxslt=1.1.39=h223e5b9_0 + - libzlib=1.3.1=hfb2fe0b_1 + - llvm-openmp=18.1.8=hde57baf_0 + - lxml=5.3.0=py39h2518886_1 + - markupsafe=2.1.5=py39h06df861_1 + - matplotlib=3.9.2=py39hdf13c20_1 + - matplotlib-base=3.9.2=py39hc57f556_1 + - matplotlib-inline=0.1.7=pyhd8ed1ab_0 + - mistune=3.0.2=pyhd8ed1ab_0 + - munkres=1.1.4=pyh9f0ad1d_0 + - nbclassic=1.1.0=pyhd8ed1ab_0 + - nbclient=0.10.0=pyhd8ed1ab_0 + - nbconvert=7.16.4=hd8ed1ab_1 + - nbconvert-core=7.16.4=pyhd8ed1ab_1 + - nbconvert-pandoc=7.16.4=hd8ed1ab_1 + - nbformat=5.10.4=pyhd8ed1ab_0 + - ncurses=6.5=hb89a1cb_0 + - nest-asyncio=1.6.0=pyhd8ed1ab_0 + - nodejs=22.9.0=h08fde81_0 + - notebook=6.4.3=pyha770c72_0 + - notebook-shim=0.2.4=pyhd8ed1ab_0 + - numpy=2.0.2=py39hd1e06cf_0 + - openjpeg=2.5.2=h9f1df11_0 + - openssl=3.3.2=h8359307_0 + - overrides=7.7.0=pyhd8ed1ab_0 + - packaging=24.1=pyhd8ed1ab_0 + - pandoc=3.2.1=hce30654_0 + - pandocfilters=1.5.0=pyhd8ed1ab_0 + - parso=0.8.4=pyhd8ed1ab_0 + - pexpect=4.9.0=pyhd8ed1ab_0 + - pickleshare=0.7.5=py_1003 + - pillow=10.4.0=py39hab9ce06_1 + - pip=24.0=pyhd8ed1ab_0 + - pkgutil-resolve-name=1.3.10=pyhd8ed1ab_1 + - platformdirs=4.2.2=pyhd8ed1ab_0 + - prometheus_client=0.20.0=pyhd8ed1ab_0 + - prompt-toolkit=3.0.47=pyha770c72_0 + - prompt_toolkit=3.0.47=hd8ed1ab_0 + - psutil=6.0.0=py39h06df861_1 + - pthread-stubs=0.4=h27ca646_1001 + - ptyprocess=0.7.0=pyhd3deb0d_0 + - pure_eval=0.2.2=pyhd8ed1ab_0 + - pycparser=2.22=pyhd8ed1ab_0 + - pygments=2.18.0=pyhd8ed1ab_0 + - pyobjc-core=10.3.1=py39hdc109a9_1 + - pyobjc-framework-cocoa=10.3.1=py39hdc109a9_1 + - pyparsing=3.1.2=pyhd8ed1ab_0 + - pysocks=1.7.1=pyha2e5f31_6 + - python=3.9.20=h9e33284_1_cpython + - python-dateutil=2.9.0=pyhd8ed1ab_0 + - python-fastjsonschema=2.20.0=pyhd8ed1ab_0 + - python-json-logger=2.0.7=pyhd8ed1ab_0 + - python_abi=3.9=5_cp39 + - pytz=2024.1=pyhd8ed1ab_0 + - pyyaml=6.0.2=py39h06df861_1 + - pyzmq=26.2.0=py39h6f9cb01_1 + - qhull=2020.2=h420ef59_5 + - qtconsole-base=5.5.2=pyha770c72_0 + - qtpy=2.4.1=pyhd8ed1ab_0 + - readline=8.2=h92ec313_1 + - referencing=0.35.1=pyhd8ed1ab_0 + - requests=2.32.3=pyhd8ed1ab_0 + - rfc3339-validator=0.1.4=pyhd8ed1ab_0 + - rfc3986-validator=0.1.1=pyh9f0ad1d_0 + - rise=5.7.1=py39h2804cbe_2 + - rpds-py=0.20.0=py39h9c3e640_1 + - scipy=1.13.1=py39h3d5391c_0 + - send2trash=1.8.3=pyh31c8845_0 + - setuptools=70.1.1=pyhd8ed1ab_0 + - six=1.16.0=pyh6c4a22f_0 + - sniffio=1.3.1=pyhd8ed1ab_0 + - soupsieve=2.5=pyhd8ed1ab_1 + - stack_data=0.6.2=pyhd8ed1ab_0 + - terminado=0.18.1=pyh31c8845_0 + - tinycss2=1.3.0=pyhd8ed1ab_0 + - tk=8.6.13=h5083fa2_1 + - tomli=2.0.1=pyhd8ed1ab_0 + - tornado=6.4.1=py39h06df861_1 + - traitlets=5.9.0=pyhd8ed1ab_0 + - types-python-dateutil=2.9.0.20240316=pyhd8ed1ab_0 + - typing-extensions=4.12.2=hd8ed1ab_0 + - typing_extensions=4.12.2=pyha770c72_0 + - typing_utils=0.1.0=pyhd8ed1ab_0 + - tzdata=2024a=h0c530f3_0 + - unicodedata2=15.1.0=py39h0f82c59_0 + - uri-template=1.3.0=pyhd8ed1ab_0 + - urllib3=2.2.2=pyhd8ed1ab_1 + - voila=0.5.7=pyhd8ed1ab_1 + - wcwidth=0.2.13=pyhd8ed1ab_0 + - webcolors=24.6.0=pyhd8ed1ab_0 + - webencodings=0.5.1=pyhd8ed1ab_2 + - websocket-client=1.8.0=pyhd8ed1ab_0 + - websockets=13.1=py39h06df861_0 + - wheel=0.43.0=pyhd8ed1ab_1 + - widgetsnbextension=4.0.11=pyhd8ed1ab_0 + - xorg-libxau=1.0.11=hb547adb_0 + - xorg-libxdmcp=1.1.3=h27ca646_0 + - xz=5.2.6=h57fd34a_0 + - yaml=0.2.5=h3422bc3_2 + - zeromq=4.3.5=hcc0f68c_4 + - zipp=3.19.2=pyhd8ed1ab_0 + - zlib=1.3.1=hfb2fe0b_1 + - zstandard=0.23.0=py39hcf1bb16_1 + - zstd=1.5.6=hb46c0d2_0 +prefix: /opt/homebrew/Caskroom/miniconda/base/envs/jupy -- GitLab