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

new requirements

parent 5aeb44a5
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 213 additions and 634 deletions
environment.yml
+ 213
− 634
View file @ a58781f3
{ name: jupy
"cells": [ channels:
{ - conda-forge
"cell_type": "code", - defaults
"execution_count": null, dependencies:
"id": "3885a60f-ddf0-4c7b-a587-6e7e037ec9b5", - anyio=4.4.0=pyhd8ed1ab_0
"metadata": { - appnope=0.1.4=pyhd8ed1ab_0
"editable": true, - argon2-cffi=23.1.0=pyhd8ed1ab_0
"slideshow": { - argon2-cffi-bindings=21.2.0=py39h06df861_5
"slide_type": "" - arrow=1.3.0=pyhd8ed1ab_0
}, - asttokens=2.4.1=pyhd8ed1ab_0
"tags": [] - async-lru=2.0.4=pyhd8ed1ab_0
}, - attrs=23.2.0=pyh71513ae_0
"outputs": [], - babel=2.14.0=pyhd8ed1ab_0
"source": [ - beautifulsoup4=4.12.3=pyha770c72_0
"import numpy as np\n", - bleach=6.1.0=pyhd8ed1ab_0
"import scipy.stats as stats\n", - brotli=1.1.0=hb547adb_1
"import matplotlib.pyplot as plt \n", - brotli-bin=1.1.0=hb547adb_1
"\n", - brotli-python=1.1.0=py39hb198ff7_1
"def f(x):\n", - bzip2=1.0.8=h93a5062_5
" return 2*x + 1\n", - ca-certificates=2024.8.30=hf0a4a13_0
"\n", - cached-property=1.5.2=hd8ed1ab_1
"n = 10\n", - cached_property=1.5.2=pyha770c72_1
"xs = np.linspace(0,4,n)\n", - certifi=2024.8.30=pyhd8ed1ab_0
"sigma_y=0.4\n", - cffi=1.17.1=py39h7f933ea_0
"ys = stats.multivariate_normal.rvs(f(xs), np.eye(n)*sigma_y**2, 1, random_state=42)\n", - charset-normalizer=3.3.2=pyhd8ed1ab_0
" \n", - comm=0.2.2=pyhd8ed1ab_0
"x_axis = np.linspace(0,4,100)\n", - contourpy=1.3.0=py39h157d57c_2
"plt.errorbar(xs,ys,yerr=sigma_y,fmt=\".\")\n", - cycler=0.12.1=pyhd8ed1ab_0
"plt.plot(x_axis, f(x_axis),'--')\n", - debugpy=1.8.6=py39hfa9831e_0
"plt.xlabel(\"x\")\n", - decorator=5.1.1=pyhd8ed1ab_0
"plt.ylabel(\"y\")\n", - defusedxml=0.7.1=pyhd8ed1ab_0
"plt.savefig(\"line.png\")\n", - entrypoints=0.4=pyhd8ed1ab_0
"plt.show()" - exceptiongroup=1.2.0=pyhd8ed1ab_2
] - executing=2.0.1=pyhd8ed1ab_0
}, - fonttools=4.54.1=py39h06df861_0
{ - fqdn=1.5.1=pyhd8ed1ab_0
"cell_type": "code", - freetype=2.12.1=hadb7bae_2
"execution_count": null, - h11=0.14.0=pyhd8ed1ab_0
"id": "9f78c037-4cca-448c-829c-dac9dea11d3c", - h2=4.1.0=pyhd8ed1ab_0
"metadata": {}, - hpack=4.0.0=pyh9f0ad1d_0
"outputs": [], - httpcore=1.0.5=pyhd8ed1ab_0
"source": [ - httpx=0.27.0=pyhd8ed1ab_0
"def chi2(x, y, sy, a, m):\n", - hyperframe=6.0.1=pyhd8ed1ab_0
" my = m * x + a\n", - icu=75.1=hfee45f7_0
" r = (y - my)/sy\n", - idna=3.7=pyhd8ed1ab_0
" return np.sum(r**2)\n", - importlib-metadata=8.0.0=pyha770c72_0
"\n", - importlib-resources=6.4.0=pyhd8ed1ab_0
"chi2(xs, ys, sigma_y, 1, 2)" - 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
"cell_type": "code", - ipython_genutils=0.2.0=pyhd8ed1ab_1
"execution_count": null, - ipywidgets=8.1.3=pyhd8ed1ab_0
"id": "76181537-7b07-4329-9615-0947d4b9dd5a", - isoduration=20.11.0=pyhd8ed1ab_0
"metadata": {}, - jedi=0.19.1=pyhd8ed1ab_0
"outputs": [], - jinja2=3.1.4=pyhd8ed1ab_0
"source": [ - json5=0.9.25=pyhd8ed1ab_0
"from matplotlib import cm\n", - jsonpointer=3.0.0=py39h2804cbe_1
"a_mesh, m_mesh = np.meshgrid(np.linspace(-3,5,100), np.linspace(-3,5,100))\n", - jsonschema=4.22.0=pyhd8ed1ab_0
"chi2_vect = np.vectorize(chi2, excluded=[0, 1])\n", - jsonschema-specifications=2023.12.1=pyhd8ed1ab_0
"m_axis = np.linspace(-3,5,100)\n", - jsonschema-with-format-nongpl=4.22.0=pyhd8ed1ab_0
"plt.plot(m_axis, chi2_vect(xs, ys, sigma_y, 1, m_axis))\n", - jupyter=1.1.1=pyhd8ed1ab_0
"plt.xlabel(\"m\")\n", - jupyter-lsp=2.2.5=pyhd8ed1ab_0
"plt.show()\n", - jupyter_client=7.4.9=pyhd8ed1ab_0
"plt.plot(m_axis, chi2_vect(xs, ys, sigma_y, m_axis, 2))\n", - jupyter_console=6.6.3=pyhd8ed1ab_0
"plt.xlabel(\"a\")\n", - jupyter_contrib_core=0.4.0=pyhd8ed1ab_0
"plt.show()\n", - jupyter_contrib_nbextensions=0.7.0=pyhd8ed1ab_0
"\n", - jupyter_core=5.7.2=pyh31011fe_1
"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", - jupyter_events=0.10.0=pyhd8ed1ab_0
"plt.show()" - 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
"cell_type": "markdown", - jupyter_server_terminals=0.5.3=pyhd8ed1ab_0
"id": "34aaaecd-6d6a-452d-8602-d8b7b9d34b72", - jupyterlab=4.2.5=pyhd8ed1ab_0
"metadata": {}, - jupyterlab_pygments=0.3.0=pyhd8ed1ab_1
"source": [ - jupyterlab_rise=0.42.0=pyhd8ed1ab_0
"## Analytisch\n", - jupyterlab_server=2.27.2=pyhd8ed1ab_0
"$\\hat m = \\frac{\\langle xy \\rangle - \\langle y \\rangle\\langle x \\rangle}{\\langle x^2 \\rangle - \\langle x \\rangle^2}$\n", - jupyterlab_widgets=3.0.11=pyhd8ed1ab_0
"\n", - kiwisolver=1.4.7=py39h157d57c_0
"$\\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", - krb5=1.21.3=h237132a_0
"\n", - lcms2=2.16=ha0e7c42_0
"\n", - lerc=4.0.0=h9a09cb3_0
"$V(\\hat m) = \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{1}{\\langle x^2 \\rangle - \\langle x \\rangle^2}$\n", - libblas=3.9.0=22_osxarm64_openblas
"\n", - libbrotlicommon=1.1.0=hb547adb_1
"\n", - libbrotlidec=1.1.0=hb547adb_1
"$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", - libbrotlienc=1.1.0=hb547adb_1
"\n", - libcblas=3.9.0=22_osxarm64_openblas
"\n", - libcxx=17.0.6=h0812c0d_3
"$\\text{cov}(\\hat m, \\hat a) = - \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{<x>}{\\langle x^2 \\rangle - \\langle x \\rangle^2}$" - libdeflate=1.20=h93a5062_0
] - libedit=3.1.20191231=hc8eb9b7_2
}, - libexpat=2.6.3=hf9b8971_0
{ - libffi=3.4.2=h3422bc3_5
"cell_type": "code", - libgfortran=5.0.0=13_2_0_hd922786_3
"execution_count": null, - libgfortran5=13.2.0=hf226fd6_3
"id": "d3b470d1-1b48-46dc-ae1c-d7a44be51e8c", - libiconv=1.17=h0d3ecfb_2
"metadata": {}, - libjpeg-turbo=3.0.0=hb547adb_1
"outputs": [], - liblapack=3.9.0=22_osxarm64_openblas
"source": [ - libopenblas=0.3.27=openmp_h517c56d_1
"C = np.eye(n)*sigma_y*sigma_y\n", - libpng=1.6.43=h091b4b1_0
"W = np.linalg.inv(C)\n", - libsodium=1.0.18=h27ca646_1
"J = np.ones((len(ys),1))\n", - libsqlite=3.46.1=hc14010f_0
"sum_xy = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*ys)))[0]\n", - libtiff=4.6.0=h07db509_3
"sum_y = (np.linalg.inv(J.T@W@J)@(J.T@W@(ys)))[0]\n", - libuv=1.48.0=h93a5062_0
"sum_x = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs)))[0]\n", - libwebp-base=1.4.0=h93a5062_0
"sum_x2 = (np.linalg.inv(J.T@W@J)@(J.T@W@(xs*xs)))[0]\n", - libxcb=1.16=hf2054a2_0
"print(sum_xy, sum_y, sum_x, sum_x2)\n", - libxml2=2.12.7=h01dff8b_4
"mhat = (sum_xy - sum_x*sum_y)/(sum_x2-sum_x*sum_x)\n", - libxslt=1.1.39=h223e5b9_0
"ahat = (sum_y*sum_x2 - sum_x*sum_xy)/(sum_x2-sum_x*sum_x)\n", - libzlib=1.3.1=hfb2fe0b_1
"print(mhat, ahat)\n", - llvm-openmp=18.1.8=hde57baf_0
"sum_siginv = sigma_y*sigma_y/len(ys)\n", - lxml=5.3.0=py39h2518886_1
"V_am = np.array([[sum_siginv/(sum_x2-sum_x*sum_x), -sum_siginv*sum_x/(sum_x2-sum_x*sum_x)],\n", - markupsafe=2.1.5=py39h06df861_1
" [-sum_siginv*sum_x/(sum_x2-sum_x*sum_x), sum_siginv*sum_x2/(sum_x2-sum_x*sum_x)]])\n", - matplotlib=3.9.2=py39hdf13c20_1
"print(V_am)\n", - matplotlib-base=3.9.2=py39hc57f556_1
"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])))" - 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
"cell_type": "markdown", - nbclient=0.10.0=pyhd8ed1ab_0
"id": "2ba98419-a39d-48e5-8ed1-47fdcab7d237", - nbconvert=7.16.4=hd8ed1ab_1
"metadata": {}, - nbconvert-core=7.16.4=pyhd8ed1ab_1
"source": [ - nbconvert-pandoc=7.16.4=hd8ed1ab_1
"# Kovarianz direkt aus C\n", - nbformat=5.10.4=pyhd8ed1ab_0
"$\\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", - ncurses=6.5=hb89a1cb_0
"\n", - nest-asyncio=1.6.0=pyhd8ed1ab_0
"$\\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" - 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
"cell_type": "code", - openjpeg=2.5.2=h9f1df11_0
"execution_count": null, - openssl=3.3.2=h8359307_0
"id": "bf423c2c-0c8f-429e-ab36-4c2d074d3575", - overrides=7.7.0=pyhd8ed1ab_0
"metadata": {}, - packaging=24.1=pyhd8ed1ab_0
"outputs": [], - pandoc=3.2.1=hce30654_0
"source": [ - pandocfilters=1.5.0=pyhd8ed1ab_0
"dmdy = sum_siginv / (sum_x2-sum_x*sum_x)*(xs-sum_x)/sigma_y**2\n", - parso=0.8.4=pyhd8ed1ab_0
"dady = sum_siginv / (sum_x2-sum_x*sum_x)*(sum_x2-xs*sum_x)/sigma_y**2\n", - pexpect=4.9.0=pyhd8ed1ab_0
"A=np.array([dmdy,dady]).T\n", - pickleshare=0.7.5=py_1003
"print(A)" - 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
"cell_type": "code", - prometheus_client=0.20.0=pyhd8ed1ab_0
"execution_count": null, - prompt-toolkit=3.0.47=pyha770c72_0
"id": "58615ab8-b1a8-433e-9380-666185592201", - prompt_toolkit=3.0.47=hd8ed1ab_0
"metadata": {}, - psutil=6.0.0=py39h06df861_1
"outputs": [], - pthread-stubs=0.4=h27ca646_1001
"source": [ - ptyprocess=0.7.0=pyhd3deb0d_0
"C_am = A.T@C@A\n", - pure_eval=0.2.2=pyhd8ed1ab_0
"print(C_am)" - 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
"cell_type": "code", - pyparsing=3.1.2=pyhd8ed1ab_0
"execution_count": null, - pysocks=1.7.1=pyha2e5f31_6
"id": "fe255a0f-c44d-4fc0-ab4b-e4bcb7c16f70", - python=3.9.20=h9e33284_1_cpython
"metadata": {}, - python-dateutil=2.9.0=pyhd8ed1ab_0
"outputs": [], - python-fastjsonschema=2.20.0=pyhd8ed1ab_0
"source": [ - python-json-logger=2.0.7=pyhd8ed1ab_0
"x_axis = np.linspace(0,4,100)\n", - python_abi=3.9=5_cp39
"plt.errorbar(xs,ys,yerr=sigma_y,fmt=\".\")\n", - pytz=2024.1=pyhd8ed1ab_0
"plt.plot(x_axis, f(x_axis),'--')\n", - pyyaml=6.0.2=py39h06df861_1
"plt.plot(x_axis, x_axis*mhat + ahat,'-')\n", - pyzmq=26.2.0=py39h6f9cb01_1
"plt.xlabel(\"x\")\n", - qhull=2020.2=h420ef59_5
"plt.ylabel(\"y\")\n", - qtconsole-base=5.5.2=pyha770c72_0
"plt.savefig(\"line.png\")\n", - qtpy=2.4.1=pyhd8ed1ab_0
"plt.show()" - readline=8.2=h92ec313_1
] - referencing=0.35.1=pyhd8ed1ab_0
}, - requests=2.32.3=pyhd8ed1ab_0
{ - rfc3339-validator=0.1.4=pyhd8ed1ab_0
"cell_type": "markdown", - rfc3986-validator=0.1.1=pyh9f0ad1d_0
"id": "b2b90d30-4e37-4ecc-866d-a45b2243d8bc", - rise=5.7.1=py39h2804cbe_2
"metadata": {}, - rpds-py=0.20.0=py39h9c3e640_1
"source": [ - scipy=1.13.1=py39h3d5391c_0
"# Fehler auf Ausgleichsgerade\n", - send2trash=1.8.3=pyh31c8845_0
"\n", - setuptools=70.1.1=pyhd8ed1ab_0
"$y = \\hat m x +\\hat a$ \n", - six=1.16.0=pyh6c4a22f_0
"\n", - sniffio=1.3.1=pyhd8ed1ab_0
"$\\frac{dy}{dm} = x$ und $\\frac{dy}{da} = 1$" - soupsieve=2.5=pyhd8ed1ab_1
] - stack_data=0.6.2=pyhd8ed1ab_0
}, - terminado=0.18.1=pyh31c8845_0
{ - tinycss2=1.3.0=pyhd8ed1ab_0
"cell_type": "code", - tk=8.6.13=h5083fa2_1
"execution_count": null, - tomli=2.0.1=pyhd8ed1ab_0
"id": "2d61e62e-77c1-41ca-a0e9-2f58bb44de08", - tornado=6.4.1=py39h06df861_1
"metadata": {}, - traitlets=5.9.0=pyhd8ed1ab_0
"outputs": [], - types-python-dateutil=2.9.0.20240316=pyhd8ed1ab_0
"source": [ - typing-extensions=4.12.2=hd8ed1ab_0
"def err(x, V):\n", - typing_extensions=4.12.2=pyha770c72_0
" A = np.array([x, 1])\n", - typing_utils=0.1.0=pyhd8ed1ab_0
" return np.sqrt(A.T@V@A)\n", - tzdata=2024a=h0c530f3_0
"\n", - unicodedata2=15.1.0=py39h0f82c59_0
"err(0, V_am)" - 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
"cell_type": "code", - webcolors=24.6.0=pyhd8ed1ab_0
"execution_count": null, - webencodings=0.5.1=pyhd8ed1ab_2
"id": "c2f71b37-b37c-4802-a356-e57d9c8df0de", - websocket-client=1.8.0=pyhd8ed1ab_0
"metadata": {}, - websockets=13.1=py39h06df861_0
"outputs": [], - wheel=0.43.0=pyhd8ed1ab_1
"source": [ - widgetsnbextension=4.0.11=pyhd8ed1ab_0
"x_axis = np.linspace(0,4,100)\n", - xorg-libxau=1.0.11=hb547adb_0
"plt.errorbar(xs,ys,yerr=sigma_y,fmt=\".\")\n", - xorg-libxdmcp=1.1.3=h27ca646_0
"plt.plot(x_axis, f(x_axis),'--')\n", - xz=5.2.6=h57fd34a_0
"plt.plot(x_axis, x_axis*mhat + ahat,'-g')\n", - yaml=0.2.5=h3422bc3_2
"plt.plot(x_axis, x_axis*mhat + ahat - np.vectorize(err,excluded=[1])(x_axis, V_am),'--g')\n", - zeromq=4.3.5=hcc0f68c_4
"plt.plot(x_axis, x_axis*mhat + ahat + np.vectorize(err,excluded=[1])(x_axis, V_am),'--g')\n", - zipp=3.19.2=pyhd8ed1ab_0
"\n", - zlib=1.3.1=hfb2fe0b_1
"plt.xlabel(\"x\")\n", - zstandard=0.23.0=py39hcf1bb16_1
"plt.ylabel(\"y\")\n", - zstd=1.5.6=hb46c0d2_0
"plt.savefig(\"line.png\")\n", prefix: /opt/homebrew/Caskroom/miniconda/base/envs/jupy
"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
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment