From da861f4f98209ae37b22ba1d845e4d039670c861 Mon Sep 17 00:00:00 2001
From: Mirko Birbaumer <mirko.birbaumer@hslu.ch>
Date: Wed, 5 Mar 2025 08:34:30 +0000
Subject: [PATCH] =?UTF-8?q?Name=20von=20Verzeichnis=20Beta=20Verteilung=20?=
=?UTF-8?q?ge=C3=A4ndert?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
notebooks/Beta Distribution/BT_2_2.ipynb | 389 -------------------
notebooks/Beta Distribution/BV_4_3.ipynb | 57 ---
notebooks/Beta Distribution/BV_4_6.ipynb | 85 ----
notebooks/Beta Distribution/beta_commands.py | 156 --------
4 files changed, 687 deletions(-)
delete mode 100644 notebooks/Beta Distribution/BT_2_2.ipynb
delete mode 100644 notebooks/Beta Distribution/BV_4_3.ipynb
delete mode 100644 notebooks/Beta Distribution/BV_4_6.ipynb
delete mode 100644 notebooks/Beta Distribution/beta_commands.py
diff --git a/notebooks/Beta Distribution/BT_2_2.ipynb b/notebooks/Beta Distribution/BT_2_2.ipynb
deleted file mode 100644
index f065766..0000000
--- a/notebooks/Beta Distribution/BT_2_2.ipynb
+++ /dev/null
@@ -1,389 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 8,
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "%matplotlib inline\n",
- "import matplotlib.pyplot as plt\n",
- "import numpy as np\n",
- "\n",
- "# Nur zur besseren Darstellung\n",
- "plt.rcParams.update({\n",
- " \"text.usetex\": True,\n",
- " \"font.family\": \"sans-serif\",\n",
- " \"font.sans-serif\": [\"Helvetica\"]})\n",
- "## for Palatino and other serif fonts use:\n",
- "plt.rcParams.update({\n",
- " \"text.usetex\": True,\n",
- " \"font.family\": \"serif\",\n",
- " \"font.serif\": [\"Palatino\"],\n",
- " \"font.size\": 8\n",
- "})"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Einleitung: Prior-Posterior\n",
- "\n",
- "\n",
- "Bei einer Münze sei die Wahrscheinlichkeit, Kopf $ K $ zu werfen, $ \\theta $. Die Wahrscheinlichkeit, Zahl $ Z $ zu werfen, ist dann $ 1-\\theta $. \n",
- "\n",
- "Nun betrachten wir einen neuen Einfränkler, der gerade von der Nationalbank in Umlauf gebracht wurde. Die Frage ist nun, wie gross ist $ \\theta $? Weil die Münze von der seriösen Nationalbank kommt, neu ist und symmetrisch aussieht, d\\\"urfen wir annehmen, dass die Münze _fair_ ist, also $ \\theta=0.5 $. \n",
- "\n",
- "Nun ist eine _reale_ Münze nie fair in dem Sinne, dass $ \\theta $ _exakt_ $ 0.5 $ ist. Es gilt vielleicht \n",
- "\n",
- "$$\n",
- "0.48<\\theta <0.52\n",
- "$$\n",
- "\n",
- "aber nie $ \\theta=0.500000000\\ldots $. Man beachte, dass es unendlich viele Werte für $ \\theta $ gibt, die zwischen $ 0 $ und $ 1 $ liegen. \n",
- "\n",
- "Wir möchten eine Aussage aufgrund von Daten machen, wie gross $ \\theta $ tatsächlich ist. Dazu werfen wir die Münze vorläufig einmal und erhalten $ K $, was wir mit $ y=1 $ bezeichnen (wir werden dies später verallgemeinern). Das heisst, wir haben Daten über die Münze gesammelt und wollen wissen, was wir über $ \\theta $ aussagen können. Wir suchen also die bedingte Wahrscheinlichkeit \n",
- "\n",
- "$$\n",
- "P(\\theta|y=1)\n",
- "$$\n",
- "\n",
- "und können die Formel von Bayes anwenden:\n",
- "\n",
- "$$\n",
- "P(\\theta|y=1)=\\frac{P(y=1|\\theta)\\cdot P(\\theta)}{P(y=1)}\n",
- "$$\n",
- "\n",
- "Nun haben wir auf der rechten Seite drei Grössen, die unbekannt sind:\n",
- "\n",
- "$$\n",
- "P(y=1|\\theta),\\qquad P(\\theta) \\quad \\text{und}\\quad P(y=1)\n",
- "$$\n",
- "\n",
- "Obwohl $ \\theta $ unendlich viele Werte annehmen kann, wollen wir $ \\theta $ vorübergehend diskretisieren. Wir nehmen an, dass $ \\theta $ nur die Werte von $ 0,0.1,0.2,\\ldots,0.9,1 $ annehmen kann. Dann gilt für die Likelihood-Funktion $ P(y=1 | \\theta) $, dass sie abhängig von $ \\theta $ ist:\n",
- "\n",
- "$$\n",
- "P(y=1|\\theta)\n",
- "=\\theta\n",
- "$$\n",
- "\n",
- "das heisst, wenn die Wahrscheinlichkeit, $ K $ zu werfen, $ \\theta $ beträgt (erstes $ \\theta $), so erscheint mit einer Wahrscheinlichkeit von $ \\theta $ (zweites $ \\theta $) auch $ K $. \n",
- "\n",
- "Ist zum Beispiel $ \\theta=0.8 $, so ist\n",
- "\n",
- "$$\n",
- "P(y=1|0.8)\n",
- "=0.8\n",
- "$$\n",
- "\n",
- "\n",
- "die Wahrscheinlichkeit, bei einer Wurfwahrscheinlichkeit von $ 0.8 $ Kopf zu werfen, natürlich auch $ 0.8 $ ist.\n",
- "\n",
- "Da wir mehrere $ \\theta $'s haben, haben wir auch mehrere $ P(y=1 | \\theta) $'s"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 9,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "[0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n",
- "[0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n"
- ]
- }
- ],
- "source": [
- "theta = np.arange(11)/10\n",
- "print(theta)\n",
- "likeli = np.linspace(0,1,num=11)\n",
- "print(likeli)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 10,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "data": {
- "text/plain": [
- "Text(50.88522208385221, 0.5, '$P(y=1\\\\mid\\\\theta )$')"
- ]
- },
- "execution_count": 10,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "plt.bar(theta,likeli, width=.01)\n",
- "plt.title(\"Likelihood\")\n",
- "plt.xlabel(r\"$\\theta$\")\n",
- "plt.ylabel(r\"$P(y=1\\mid\\theta )$\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Nun kommen wir zur Prior-Verteilung $ P(\\theta) $, und hier gibt es sehr viele Möglichkeiten. Diese müssen wir allerdings vorher festlegen. Denken wir uns die Münze als eher fair, wissen es aber nicht genau, können wir eine Wahrscheinlichkeitsverteilung wie folgt w\\\"ahlen.\n",
- "\n",
- "Die Prior-Werte für die $ \\theta $'s sind"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 11,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "[0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n",
- "[0 1 2 3 4 5 4 3 2 1 0]\n",
- "[0. 0.04 0.08 0.12 0.16 0.2 0.16 0.12 0.08 0.04 0. ]\n"
- ]
- }
- ],
- "source": [
- "theta = np.arange(11)/10\n",
- "y = np.array([0,1,2,3,4,5,4,3,2,1,0])\n",
- "prior = y / np.sum(y)\n",
- "print(theta)\n",
- "print(y)\n",
- "print(prior)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 12,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "data": {
- "text/plain": [
- "Text(50.88522208385221, 0.5, '$P(\\\\theta)$')"
- ]
- },
- "execution_count": 12,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "plt.bar(theta,prior, width=.01)\n",
- "plt.title(\"Prior\")\n",
- "plt.xlabel(r\"$\\theta$\")\n",
- "plt.ylabel(r\"$P(\\theta)$\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Nun können wir die Posterior-Wahrscheinlichkeiten für alle $ \\theta $'s, also die Posterior-Verteilung, berechnen. So ist zum Beispiel\n",
- "\n",
- "$$\n",
- "P(0.8 | y=1)\n",
- "=\\frac{P(y=1 | 0.8)\\cdot P(0.8)}{P(y=1)}\n",
- "=\\frac{0.8\\cdot 0.8}{0.5}\n",
- "=0.128\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Nun updaten wir unsere Wahrscheinlichkeiten mit der Formel von Bayes für alle $\\theta$s, mit der Beobachtung, dass Kopf geworfen wurde ($y=1$):"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 13,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Zähler [0. 0.004 0.016 0.036 0.064 0.1 0.096 0.084 0.064 0.036 0. ]\n",
- "Nenner 0.5\n",
- "Posterior [0. 0.008 0.032 0.072 0.128 0.2 0.192 0.168 0.128 0.072 0. ]\n"
- ]
- }
- ],
- "source": [
- "posterior = (likeli*prior) / np.sum(likeli*prior)\n",
- "\n",
- "print(\"Zähler\", likeli*prior)\n",
- "print(\"Nenner\", np.sum(likeli*prior))\n",
- "print(\"Posterior\", posterior)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 14,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "data": {
- "text/plain": [
- "Text(50.88522208385221, 0.5, '$P(\\\\theta\\\\mid y=1)$')"
- ]
- },
- "execution_count": 14,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "plt.bar(theta,posterior, width=.01)\n",
- "plt.title(\"Posterior\")\n",
- "plt.xlabel(r\"$\\theta$\")\n",
- "plt.ylabel(r\"$P(\\theta\\mid y=1)$\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Wir sehen, dass die Prior-Verteilung durch eine leicht andere Gewichtung in die Posterior-Verteilung einfliesst. Die $ \\theta $'s über $ 0.5 $ werden stärker gewichtet als die unter $ 0.5 $. \n",
- "\n",
- "Hätten wir Zahl ($ Z $) geworfen, so hätten wir eine höhere Gewichtung für Wurfwahrscheinlichkeiten unter $ 0.5 $."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Verschiedene Priors\n",
- "\n",
- "Wir betrachten als erstes eine uniform verteilte Prior-Verteilung."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 15,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "[0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n",
- "[0 1 1 1 1 1 1 1 1 1 0]\n",
- "[0. 0.11111111 0.11111111 0.11111111 0.11111111 0.11111111\n",
- " 0.11111111 0.11111111 0.11111111 0.11111111 0. ]\n"
- ]
- },
- {
- "data": {
- "text/plain": [
- "Text(50.88522208385221, 0.5, '$P(\\\\theta)$')"
- ]
- },
- "execution_count": 15,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "theta = np.arange(11)/10\n",
- "y = np.array([0,1,1,1,1,1,1,1,1,1,0])\n",
- "prior = y / np.sum(y)\n",
- "print(theta)\n",
- "print(y)\n",
- "print(prior)\n",
- "plt.bar(theta,prior, width=.01)\n",
- "plt.title(\"Prior\")\n",
- "plt.xlabel(r\"$\\theta$\")\n",
- "plt.ylabel(r\"$P(\\theta)$\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 16,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "[0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n",
- "[ 0. 1. 1.2 2. 5. 15. 5. 2.25 1.2 1. 0. ]\n",
- "[0. 0.02971768 0.03566122 0.05943536 0.14858841 0.44576523\n",
- " 0.14858841 0.06686478 0.03566122 0.02971768 0. ]\n"
- ]
- },
- {
- "data": {
- "text/plain": [
- "Text(50.88522208385221, 0.5, '$P(\\\\theta)$')"
- ]
- },
- "execution_count": 16,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "theta = np.arange(11)/10\n",
- "y = np.array([0,1,1.2,2,5,15,5,2.25,1.2,1,0])\n",
- "prior = y / np.sum(y)\n",
- "print(theta)\n",
- "print(y)\n",
- "print(prior)\n",
- "plt.bar(theta,prior, width=.01)\n",
- "plt.title(\"Prior\")\n",
- "plt.xlabel(r\"$\\theta$\")\n",
- "plt.ylabel(r\"$P(\\theta)$\")"
- ]
- }
- ],
- "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.10.11"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
diff --git a/notebooks/Beta Distribution/BV_4_3.ipynb b/notebooks/Beta Distribution/BV_4_3.ipynb
deleted file mode 100644
index 5d75b09..0000000
--- a/notebooks/Beta Distribution/BV_4_3.ipynb
+++ /dev/null
@@ -1,57 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 13,
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "%matplotlib inline\n",
- "import beta_commands as bt\n",
- "import matplotlib.pyplot as plt\n",
- "plt.rcParams['figure.figsize'] = [12,14]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 17,
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "bt.prior(a=5.8, b=4.2, z=4, N=10) "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "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.10.11"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
diff --git a/notebooks/Beta Distribution/BV_4_6.ipynb b/notebooks/Beta Distribution/BV_4_6.ipynb
deleted file mode 100644
index e2d130c..0000000
--- a/notebooks/Beta Distribution/BV_4_6.ipynb
+++ /dev/null
@@ -1,85 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Bestimmung 95%-HDI für die Betaverteilung\n",
- "\n",
- "Die folgende Prozedur definiert, das HDI für die Beta-Verteilung:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 4,
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "from scipy.stats import beta\n",
- "import numpy as np\n",
- "\n",
- "def hdi(a,b, prob = 0.95):\n",
- " k = 0\n",
- " x = np.linspace(0,1,1000)\n",
- " y = beta.pdf(x,a,b)\n",
- " while True:\n",
- " k = k+0.0001\n",
- " if np.sum(y[y > k])/np.size(x) < prob:\n",
- " break\n",
- " return x[np.argwhere(y > k)][0] ,x[np.argwhere(y > k)][np.argwhere(y > k).size-1] "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "In `hdi(...)` müssen nur noch die Parameter $a$ und $b$ eingegeben werden. Wenn die Prozentzahl noch geändert soll, fügen Sie ein zusätzliches `prob=...` ein. prob=0.5 würde dann bedeuten, dass 50% der glaubwürdigsten Parameter im HDI enthalten sind. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 5,
- "metadata": {
- "tags": []
- },
- "outputs": [
- {
- "data": {
- "text/plain": [
- "(array([0.14614615]), array([0.85385385]))"
- ]
- },
- "execution_count": 5,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "hdi(a=3, b=3) "
- ]
- }
- ],
- "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.10.11"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
diff --git a/notebooks/Beta Distribution/beta_commands.py b/notebooks/Beta Distribution/beta_commands.py
deleted file mode 100644
index 0d1072b..0000000
--- a/notebooks/Beta Distribution/beta_commands.py
+++ /dev/null
@@ -1,156 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Jan 22 13:14:42 2021
-
-@author: bl
-"""
-
-from scipy.stats import beta
-
-import matplotlib.pyplot as plt
-import numpy as np
-import arviz as av
-
-
-
-
-def hdi(a,b, prob):
- k = 0
- x = np.linspace(0,1,1000)
- y = beta.pdf(x,a,b)
- while True:
- k = k+0.001
- if np.sum(y[y > k])/np.size(x) < prob:
- break
- return np.array([x[np.argwhere(y > k)][0] ,x[np.argwhere(y > k)][np.argwhere(y > k).size-1]])
-
-x = np.linspace(0,1,1000)
-
-
-
-def prior(a,b,z,N, prob=.95):
- plt.subplots(3, 1)
- plt.subplot(311)
- y = beta.pdf(x,a,b)
- omega = np.round((z+a-1)/(z+a+N-z+b-2),3)
- pr = int(prob*100)
- plt.plot(x,y,color="seagreen")
- plt.title("Prior (beta)")
- plt.xlabel("θ")
- plt.ylabel(r"beta$(θ\mid$"+str(a)+" ," +str(b)+")")
- if a>1 and b > 1:
- omega_prior = np.round((a-1)/(a+b-2),3)
- xl = np.round(hdi(a,b,prob)[0][0],3)
- xr = np.round(hdi(a,b,prob)[1][0],3)
- yl = beta.pdf(xl,a,b)
- plt.plot([xl,xl],[0,yl], linestyle="dashed",color="black")
- plt.plot([xr,xr],[0,yl], linestyle="dashed",color="black")
- plt.plot([xl,xr],[yl,yl],color="orange")
- plt.text(xl-.01,yl, str(xl),ha="right",va="bottom")
- plt.text(xr+.01,yl, str(xr),ha="left",va="bottom")
- plt.text((xr+xl)/2,yl+.1*beta.pdf(0.5,a,b),str(pr)+"% HDI",ha="center")
- plt.text(0,beta.pdf(omega,z+a,N-z+b)*.75,"Modus: "+str(omega_prior),ha="left")
-
- elif a==1 and b==1:
- plt.text(0,beta.pdf(omega,z+a,N-z+b)*.75,"Modus: 0.5",ha="left")
- plt.ylim(0,beta.pdf(omega,z+a,N-z+b))
- plt.fill_between(x,y,facecolor="mediumaquamarine")
-
-
- plt.subplot(312)
-
- y = x**z*(1-x)**(N-z)
- plt.ylim(0, np.max(y)*1.1)
- plt.plot(x,y,color="seagreen")
- plt.title("Likelihood (Bernoulli)")
- plt.xlabel(r"$θ$")
- plt.ylabel(r"$p(D\mid θ)$")
- plt.fill_between(x,y,facecolor="mediumaquamarine")
- plt.text(0, np.max(y)*.75, "Daten: $z=\ $"+str(z)+", $N =\ $"+str(N),ha="left")
- plt.text(0, np.max(y)*.4, "Max bei 0.85",ha="left")
-
-
-
-
-
- plt.subplot(313)
- a, b = z+a, N-z+b
- y = beta.pdf(x,a,b)
- plt.plot(x,y,color="seagreen")
- plt.title("Posterior (beta)")
- plt.xlabel(r"$θ$")
- plt.ylabel(r"beta$(θ\mid$"+str(a)+" ," +str(b)+")")
- plt.fill_between(x,y,facecolor="mediumaquamarine")
- xl = np.round(hdi(a,b,prob)[0][0],3)
- xr = np.round(hdi(a,b,prob)[1][0],3)
- yl = beta.pdf(xl,a,b)
- plt.plot([xl,xl],[0,yl], linestyle="dashed",color="black")
- plt.plot([xr,xr],[0,yl], linestyle="dashed",color="black")
- plt.plot([xl,xr],[yl,yl],color="orange")
- plt.text(xl-.01,yl, str(xl),ha="right",va="bottom")
- plt.text(xr+.01,yl, str(xr),ha="left",va="bottom")
-
- plt.text((xr+xl)/2,yl+.1*beta.pdf(omega, a,b),str(pr)+"% HDI",ha="center")
-
- plt.text(0,beta.pdf(omega, a,b)*.75,"Modus: "+str(omega),ha="left")
- plt.fill_between(x,y,facecolor="mediumaquamarine")
- plt.ylim(0,beta.pdf(omega, a,b))
- plt.tight_layout(w_pad=3, h_pad=3)
-
-
-def beta_hist(a,b,N):
- x = np.linspace(0,1,1000)
-
- plt.subplots(1, 2, sharey=True,figsize=(10,5))
-
- plt.subplot(121)
- y = beta.pdf(x,a,b)
- plt.plot(x,y,color="seagreen")
- plt.title("Exakte Verteilung")
- plt.xlabel("θ")
- plt.ylabel("p(θ)")
- xl = np.round(hdi(a,b)[0][0],3)
- xr = np.round(hdi(a,b)[1][0],3)
-
- yl = beta.pdf(xl,a,b)
- plt.plot([xl,xl],[0,yl], linestyle="dashed",color="black")
- plt.plot([xr,xr],[0,yl], linestyle="dashed",color="black")
- plt.plot([xl,xr],[yl,yl],color="orange")
- plt.text(xl-.01,yl, str(xl),ha="right",va="bottom")
- plt.text(xr+.01,yl, str(xr),ha="left",va="bottom")
- omega = np.round((a-1)/(a+b-2),3)
- plt.text(0,3.5,"Modus: "+str(omega),ha="left")
- plt.text((xr+xl)/2,yl+.075*beta.pdf(omega, a,b),r"95$\%$ HDI",ha="center")
- plt.fill_between(x,y,facecolor="mediumaquamarine")
- plt.ylim(0,4)
-
- plt.subplot(122)
-
- y = beta.rvs(a,b,size=N)
- bin = np.linspace(0,1,51)
-
- plt.title("N="+str(N))
- plt.xlabel("θ")
- plt.ylabel("p(θ)")
- n, b = np.histogram(y,bins=bin,density=True)
- bin_max = np.where(n == n.max())
- omega = np.round(b[bin_max][0],2)
- xl = np.round(av.hdi(y,hdi_prob=0.95)[0],3)
- xr = np.round(av.hdi(y,hdi_prob=0.95)[1],3)
- yl = .5
- plt.hist(y,bins=bin,edgecolor="black",density=True, color="mediumaquamarine",zorder=1)
- plt.plot([xl,xl],[0,yl], linestyle="dashed",color="black")
- plt.plot([xr,xr],[0,yl], linestyle="dashed",color="black")
- plt.plot([xl,xr],[yl,yl],color="orange")
- plt.text(xl-.01,yl, str(xl),ha="right",va="bottom")
- plt.text(xr+.01,yl, str(xr),ha="left",va="bottom")
- plt.text(0,3.5,"Modus: "+str(omega),ha="left")
- plt.text((xr+xl)/2,yl+.075,"95% HDI",ha="center")
- plt.ylim(0,4)
-
-
-
-
-
-
--
GitLab