pgp1-python-einfuehrung/Herleitung_Methode_der_kleinsten_Quadrate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f3b16f1a-521e-4caf-828e-85251faf2c6c",
   "metadata": {},
   "source": [
    "Dieses Dokument beinhaltet die Herleitung der Methode der kleinsten Quadrate, die im Kapitel 1 beim Fitten benutzt wird."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a144903b-e362-40de-8af0-7da7b39fe260",
   "metadata": {},
   "source": [
    "### Methode der kleinsten Quadrate\n",
    "\n",
    "Im folgenden wolllen wir die **Methode der kleinsten Quadrate (Least Squares)** näher beleuchten. Diese Methode wird oft benutzt, um eine Funktion $\\lambda(x; \\ $**$\\phi$**$)$ mit den Funktionsparametern $\\mathbf{\\phi}$ an die gemessenen Punkte **$(x,y)$** anzupassen. Um jedoch die **Methode der kleinsten Quadrate** zu verstehen, wollen wir sie erst einmal anschaulich und mathematisch herleiten. Dabei stüzen wir uns im Folgenden auf eine Herleitung aus dem Buch **\"Statistical Data Analysis\"**  von **Glen Cowan**.\n",
    "\n",
    "In unserem Grundpraktikum haben wir bereits gelernt, dass Messwerte durch Zufallszahlen $x_i$ representiert werden und einer gewissen **Wahrscheinlichkeitsdichtefunktion (probability density function)** $f(x)$ unterliegen.\n",
    "\n",
    "<figure class=\"image\">\n",
    "<img src=\"images/MaterialPythonkurs092018/PorbDensFun.png\"  alt=\"{{ Beispiel PDF }}\" width=70%>\n",
    "</figure>\n",
    "\n",
    "\n",
    "Eine **pdf** gibt an, mit welcher **Wahrscheinlichkeit ein Wert $x_i$** innerhalb eines **infinitesimalen Intervals $\\text{d}x_i$** zu finden ist. Des Weitren gilt, dass die Gesamtwahrscheinlichkeit gegeben ist durch $\\int_S f(x) dx = 1$. \n",
    "\n",
    "Nun betrachten wir folgendes Beispiel: In unserem Labor messen wir genau drei mal die Raumtemperartur T. Auch hier gilt, dass unsere Messung der einzelnen $T_i$ einer gewissen **Wahrscheinlichkeitsdichtefunktion** folgen. Betrachten Sie nun das folgende Bild; Welche **Wahrscheinlichkeitsdichtefunktion** passt besser zu den gezeigten Daten und **Warum?**\n",
    "\n",
    "<figure class=\"image\">\n",
    "<img src=\"images/MaterialPythonkurs092018/ProbMaxTemp.png\"  alt=\"{{ Beispiel PDF }}\" width=100%>\n",
    "</figure>\n",
    "\n",
    "Die rechte Verteilung spiegelt unsere Messdaten besser wieder. Dies können wir auch mathematisch ausdrücken. Für $N$ voreinander unabhängige Zufallszahlen bzw. Messpunkte (in unserem Beispiel $N = 3$) ist die Gesamtwahrscheinlichkeit gegeben durch das Produkt der einzelnen Wahrscheinlichkeitsdichten $f(x_i, \\theta)$ multipliziert mit dem jeweiligen infinitesimalen element $dx_i$\n",
    "\n",
    "$$\\prod_{i = 1}^{N} f(x_i,\\theta) \\  dx_i \\text{   für alle } x_i \\text{ in } [x_i, x_i + dx_i]$$\n",
    "\n",
    "wobei $x_i$ in unserem Beispiel den Messpunkten $T_i$ und $f(x_i,\\theta)$ unserer Gausverteilung mit $\\theta = (\\mu, \\sigma)$ entspricht. Sprich sofern unsere Werte gut von der jeweiligen **Wahrscheinlichkeitsdichtefunktion** repräsentiert werden, d.h. wir die richtigen Parameter $\\theta$ gewählt haben (wie im rechten oberen Plot), gilt \n",
    "\n",
    "$$ \\prod_{i = 1}^{N} f(x_i,\\theta)  dx_i$$ \n",
    "\n",
    "ist **maximal**. Da die einzelnen $dx_i$ von unseren Parametern $\\theta$ unabhängig sind, gilt die gleiche Argumentation auch für \n",
    "\n",
    "$$ \\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N) = \\prod_{i = 1}^{N} f(x_i,\\theta)$$ \n",
    "\n",
    "wobei $\\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N)$ die sogenannte **likely hood function** darstellt.\n",
    "\n",
    "Wie kommen wir nun von der **likely hood function** auf unsere **Methode der kleinsten Quadrate** und das Fitten einer Funktion $\\lambda(x; \\ $**$\\phi$**$)$ an die gemessenen Punkte **$(x,y)$**? Dazu brauche wir noch einen Zwischenschritt. Oftmals ist es einfacher, statt die **likely hood function** zu maximieren, die so genannte **log likely hood function**\n",
    "\n",
    "$$ \\log( \\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N)) = \\sum_{i = 1}^{N} \\log(f(x_i,\\theta))$$\n",
    "\n",
    "zu maximieren. Dies ist im Grunde das Gleiche, da der Logarithmus eine monoton-steigende Funktion ist. Auch in unserem Fall der **Methode der kleinsten Quadrate** benötigen wir die **log likely hood function**. \n",
    "\n",
    "Stellen Sie sich nun vor, wir haben eine Messung mit $N$ voneinander unabhängigen Messpunkten (x,y). Des Weiteren nehmen wir an, dass alle $x_i$ ohne Fehler sind und dass unsere $y_i$ gaußförmig um einen unbekannten Wahrenwert $\\lambda_i$ (sprich $\\lambda_i$ entspricht dem Erwartungswert $\\mu_i$ unserer Gaußverteilung) mit einer bekannten Varianz $\\Delta y_i^2$ verteilt sind (Diese Annahme lässt sich mit dem zentralen Grenzwertsatz begründen, so lange der Fehler sich aus der Summe kleiner Fehler zusammensetzt). Die dazugehörige  **likely hood function** ist dann gegeben durch:\n",
    "\n",
    "$$ \\mathcal{L}(y_1 ... y_N; \\lambda_1 ... \\lambda_N, \\Delta y_1 ... \\Delta y_N)) = \\prod_{i = 1}^{N}\\frac{1}{\\sqrt{2 \\pi \\Delta y_i^2}} \\cdot \\exp \\bigg( \\frac{ -(y_i - \\lambda_i)^2}{2 \\cdot \\Delta y_i^2}\\bigg)$$\n",
    "\n",
    "Beziehungsweise die **log likely hood function** mit $\\lambda_i = \\lambda(x_i; \\phi)$ ergibt sich zu\n",
    "\n",
    "$$ \\log(\\mathcal{L}(y, \\theta)) \\approx -\\frac{1}{2} \\sum_{i = 1}^{N}\\bigg( \\frac{ (y_i - \\lambda(x_i; \\phi))^2}{\\Delta y_i^2}\\bigg)$$\n",
    "\n",
    "wobei die konstanten Terme, welche nicht von unserer Funktion $\\lambda(x_i; \\phi)$ abhängen, vernachlässigt worden sind. Durch den Faktor $-\\frac{1}{2}$ ist das Maximieren dieser **log likely hood function** gleich dem Minimieren von\n",
    "\n",
    "$$ \\chi(\\phi_1 ... \\phi_N)^2 = \\sum_{i = 1}^{N} \\frac{ (y_i - \\lambda(x_i; \\phi))^2}{\\Delta y_i^2}$$\n",
    "\n",
    "Diese Funktion ist unsere gesuchte **Methode der kleinsten Quadrate**. Mit ihrer Hilfe kann eine beliebige Funktion $\\lambda(x; \\phi)$, welche liniear in ihren Parametern $\\phi$ ist, an unsere Messdaten $(x,y\\pm\\Delta y)$ gefittet werden. Dabei stellt der Fitprozess selbst lediglich ein Minimierungsproblem dar. Im Folgenden sind unsere Annahmen noch einmal grafisch in einem Beispiel dargestellt.\n",
    "\n",
    "<figure class=\"image\">\n",
    "<img src=\"images/MaterialPythonkurs092018/LeastSquare.png\"  alt=\"{{ Least Square Beispiel }}\" width=100%>\n",
    "</figure>\n",
    "\n",
    "Es gibt verschiedene Arten von Algorithmen um Minimierungsprobleme zu lösen. Wie diese genau aufgebaut sind, lernen Sie in anderen Progrmmierveranstaltungen wie zum Beispiel *Programmieren für Physiker* oder *Computer in der Wissenschaft*. Zum Glück haben uns bereits in Python andere Menschen diese Arbeit abgenommen und wir können aus dem Package `scipy.optimize` die Funktion `curve_fit` verwenden.\n",
    "\n",
    "Hierbei stellt `curve_fit` eine Methode dar, Fit-Funktionen nach der obigen vorgestellten Methode der *kleinsten Quadraten* zu bestimmen. Dies hat zur Folge, dass lediglich die y-Fehler der Messwerte für den Fit verwendet werden können."
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}