{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Distribuzione di Poisson\n", "========================" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La distribuzione di Poisson:\n", "$$ P_\\lambda(n) = \\frac{\\lambda^n}{n!} e^{-\\lambda}$$\n", "fornisce la probabilità di osservare $n$ eventi in un processo in cui il numero medio di osservazioni è $\\lambda$.\n", "\n", "## Proprietà:\n", "\n", "$$\\sum_{n=0}^\\infty P_\\lambda(n) = 1$$\n", "#### Valor medio:\n", "$$ = e^{-\\lambda} \\sum_{n=0}^\\infty \\frac{n \\lambda^n}{n!} = e^{-\\lambda}\\sum_{n=1}^\\infty \\frac{\\lambda^n}{(n-1)!} = \\lambda \\, e^{-\\lambda} \\sum_{n=0}^\\infty \\frac{\\lambda^n}{n!} = \\lambda $$\n", "oppure\n", "$$ = e^{-\\lambda} \\sum_{n=0}^\\infty \\frac{n \\lambda^n}{n!} = e^{-\\lambda} \\sum_{n=0}^\\infty \n", "\\frac{\\lambda}{n!} \\frac{d \\lambda^n}{d \\lambda} = \\lambda \\, e^{-\\lambda} \\frac{d }{d \\lambda}\\left(\\sum_{n=0}^\\infty \\frac{\\lambda^n}{n!} \\right)= \\lambda.$$\n", "In modo analogo si ha: $\\, \\,= \\lambda^2$.\n", "\n", "#### Scarto quadratico medio:\n", "$$ <(n-)^2> = - ^2 = e^{-\\lambda} \\sum_{n=0}^\\infty \\frac{n^2 \\lambda^n}{n!} - \\lambda^2 = \n", "e^{-\\lambda} \\sum_{n=0}^\\infty \\frac{(n(n-1) +n)\\, \\lambda^n}{n!} - \\lambda^2 = \n", "\\lambda^2 + \\lambda - \\lambda^2 = \\lambda $$\n", "oppure, essendo\n", "$$\\frac{d^2 \\lambda^n}{d \\lambda^2} = n(n-1)\\lambda^{n-2}, \\qquad \\frac{d \\lambda^n}{d \\lambda} = n\\lambda^{n-1}$$\n", "e quindi\n", "$$n^2 \\lambda^n = \\lambda^2 \\frac{d^2 \\lambda^n}{d \\lambda^2} + n \\lambda^n = \\lambda^2 \\frac{d^2 \\lambda^n}{d \\lambda^2} + \\lambda \\frac{d \\lambda^n}{d \\lambda}$$\n", "$$ e^{-\\lambda} \\sum_{n=0}^\\infty \\frac{n^2 \\lambda^n}{n!} = e^{-\\lambda} \n", "\\left( \\lambda^2 \\frac{d^2}{d \\lambda^2} \\sum_{n=0}^\\infty \\frac{\\lambda^n}{n!}\n", "+ \\lambda \\frac{d}{d \\lambda}\\sum_{n=0}^\\infty \\frac{\\lambda^n}{n!} \\right) = \\lambda^2 + \\lambda$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.special" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Definiamo una funzione che calcoli il valore della distribuzione di Poisson a fissato n e $\\lambda$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def PoissonVal(n,lam):\n", " '''\n", " Returns value of the Poisson distribution with average lam at n.\n", " Uses np.power and scipy.special.factorial to vectorize.\n", " It doesn't work for large n, the factorial becomes too big.\n", " '''\n", " return np.exp(-lam)*np.power(lam,n)/scipy.special.factorial(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quando n diventa grande, n! diventa enorme e Python (come tutti gli altri linguaggi) non riesce più a completare il calcolo.\n", "\n", "È necessario utilizzare delle formule approssimate per il fattoriale che permettono di calcolare il logaritmo della distribuzione di Poisson che viene esponenziata solo nell'ultimo passaggio.\n", "\n", "### Approssimazioni del fattoriale\n", "Approssimazione di Stirling:\n", "$$ n! \\approx \\sqrt{2\\pi n}\\,\\left(\\frac{n}{e}\\right)^n$$\n", "Approssimazione di Ramanujan:\n", "$$ n! \\approx \\sqrt{\\pi}\\,\\left(\\frac{n}{e}\\right)^n\\sqrt[6]{8 n^3 + 4 n^2 + n + \\frac{1}{30}}$$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "def BigPoissonVal(n,lam):\n", " '''\n", " Returns value of the Poisson distribution with average lam at n\n", " for large value of n using Stirling's Formula for the factorial:\n", " n! = sqrt(2*pi*n)(n/e)**n\n", " '''\n", " t=-lam+n*np.log(lam/n)+n-np.log(2*np.pi*n)/2\n", " return np.exp(t)\n", "\n", "def BigPoissonVal2(n,lam):\n", " '''\n", " Returns value of the Poisson distribution with average lam at n\n", " for large value of n using Ramanujan's Formula for the factorial\n", " '''\n", " t=-lam+n*np.log(lam/n)+n-np.log(np.pi)/2. -np.log(n*(1.+4.*n*(1.+2.*n))+1./30.)/6.\n", " return np.exp(t)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Test del valore esatto vs Stirling vs Ramanujan" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lam = 2.\n", "print(' Esatto Stirling Ramanujan')\n", "print('')\n", "for n in range(5,10):\n", " print(f\"{PoissonVal(n,lam):21.19f} {BigPoissonVal(n,lam):21.19f} {BigPoissonVal2(n,lam):21.19f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Esempi di distribuzione di Poisson" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(15)\n", "\n", "fig, ax = plt.subplots(figsize=(12, 5))\n", "ax.set_title('Poisson distribution')\n", "ax.set_xlabel('n')\n", "ax.set_ylabel('P(n)')\n", "ax.grid(axis='x')\n", "ax.plot(t,PoissonVal(t,1),marker='o',c='r', linestyle = '--',label='$\\lambda$ = 1')\n", "ax.plot(t,PoissonVal(t,4),marker='o',c='b', linestyle = '--',label='$\\lambda$ = 4')\n", "ax.plot(t,PoissonVal(t,6),marker='o',c='g', linestyle = '--',label='$\\lambda$ = 6')\n", "ax.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Altre distribuzioni interessanti collegate alla distribuzione di Poisson\n", "\n", "Definiamo due funzioni che calcolino la distribuzione normale e quella binomiale" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def GaussVal(x,mu,sigma):\n", " return np.exp(-((x-mu)/sigma)**2/2)/sigma/np.sqrt(2*np.pi)\n", "\n", "def BinomialVal(n,ntot,p):\n", " '''\n", " Returns value of the Binomial distribution for n successes in ntot tries with\n", " probability p for a success.\n", " Uses np.power and scipy.special.factorial to vectorize\n", " '''\n", " n1 = ntot-n\n", " factor = (scipy.special.factorial(ntot)/\n", " scipy.special.factorial(n)/scipy.special.factorial(n1))\n", " return np.power(p,n)*np.power(1-p,n1)*factor" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Descrizione del decadimento di nuclei radioattivi usando la distribuzione di Poisson e la distribuzione binomiale\n", "\n", "Supponiamo di avere N nuclei con probabilità P di decadere nell'unità di tempo T0. Il valor medio dei conteggi in T0 sarà ave = N\\*P
\n", "Generiamo ntry pseudoesperimenti: creiamo ntry array di numeri casuali x fra 0 e 1.25 di dimensione N. I valori per cui x > 1.0 individuano i nuclei che decadono." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nuclei = 150\n", "ntry = 50000\n", "prob_decay = 0.05\n", "# Ho 1-p eventi rifiutati nell'intervallo [0,1]. \n", "# Per selezionare eventi con probabilita' p scegliendo quelli maggiori di uno\n", "# devo scalare l'intervallo per alpha=1/(1-p)\n", "# E.g: p=0.2 -> alpha = 1/0.8 = 1.25\n", "PseudoExp = np.random.uniform(size=(ntry,nuclei))/(1.-prob_decay)\n", "#PseudoExp\n", "\n", "#Check PseudoExp\n", "nbins = 25\n", "xrange = (0.,1./(1.-prob_decay))\n", "figc, axc = plt.subplots()\n", "nevent, bins, patches = axc.hist(np.ndarray.flatten(PseudoExp), nbins, range=xrange)\n", "print(bins)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Select values larger than 1 by using the floor method\n", "PseudoExp1 = np.floor(PseudoExp)\n", "print(PseudoExp1)\n", "# Count the decay taking the sum over lines\n", "decays = np.sum(PseudoExp1,axis=1).astype(int)\n", "decays" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Select values larger than 1\n", "# Gives the same result of the previous method\n", "#for i in range(ntry):\n", "# temp = PseudoExp[i]\n", "# temp_decayed = [x for x in temp if x > 1]\n", "# print(len(temp_decayed))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ave = decays.mean()\n", "variance = decays.var()\n", "std_dev = decays.std()\n", "print(\"average counts:\",ave)\n", "print(\"standard deviation:\",std_dev)\n", "print(\"standard deviation squared:\",std_dev**2)\n", "print(\"variance:\",variance)\n", "unique, counts = np.unique(decays, return_counts=True)\n", "print(\"unique values:\",unique)\n", "print(\"counts:\",counts)\n", "total = np.sum(counts)\n", "counts = counts/total\n", "print(\"counts:\",counts)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(3*np.floor(ave))\n", "average_counts = nuclei*prob_decay\n", "yP = PoissonVal(t,average_counts)\n", "yB = BinomialVal(t,nuclei,prob_decay)\n", "\n", "# check normalization \n", "print('Poisson Normalization:',np.sum(yP))\n", "\n", "fig1, ax1 = plt.subplots(figsize=(12, 5))\n", "ax1.set_title('Poisson distribution vs data for $\\lambda$ ='+str(average_counts))\n", "ax1.set_xlabel('n')\n", "ax1.set_ylabel('P(n)')\n", "ax1.grid(axis='x')\n", "#leg1 = '$\\lambda$ ='+str(average_counts)\n", "ax1.plot(t,yP,marker='o',c='b', linestyle = '--',label='$\\lambda$ ='+str(average_counts))\n", "ax1.plot(t,yB,marker='o',c='k', linestyle = '--',label='Binomial(n,'+str(nuclei)+','+str(prob_decay)+')')\n", "#ax1.scatter(unique,counts,marker='o',c='g',label='data')\n", "ax1.plot(unique,counts,marker='o',c='g',label='data')\n", "ax1.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Caso complicato con prob_decay = 0.2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nuclei = 150\n", "ntry = 10000\n", "prob_decay = 0.2\n", "PseudoExp = np.random.uniform(size=(ntry,nuclei))/(1.-prob_decay)\n", "\n", "# Select values larger than 1 by using the floor method\n", "PseudoExp1 = np.floor(PseudoExp)\n", "print(PseudoExp1)\n", "# Count the decay taking the sum over lines\n", "decays = np.sum(PseudoExp1,axis=1).astype(int)\n", "\n", "ave = decays.mean()\n", "std_dev = decays.std()\n", "print(\"average counts:\",ave)\n", "print(\"standard deviation:\",std_dev)\n", "print(\"standard deviation squared:\",std_dev**2)\n", "unique, counts = np.unique(decays, return_counts=True)\n", "total = np.sum(counts)\n", "counts = counts/total" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(3*np.floor(ave))\n", "average_counts = nuclei*prob_decay\n", "yP = PoissonVal(t,average_counts)\n", "yB = BinomialVal(t,nuclei,prob_decay)\n", "# check normalization \n", "print('Poisson Normalization:',np.sum(yP))\n", "\n", "fig1, ax1 = plt.subplots(figsize=(12, 5))\n", "ax1.set_title('Poisson distribution vs data for $\\lambda$ ='+str(average_counts))\n", "ax1.set_xlabel('n')\n", "ax1.set_ylabel('P(n)')\n", "ax1.grid(axis='x')\n", "#leg1 = '$\\lambda$ ='+str(average_counts)\n", "ax1.plot(t,yP,marker='o',c='b', linestyle = '--',label='$\\lambda$ ='+str(average_counts))\n", "ax1.plot(t,yB,marker='o',c='k', linestyle = '--',label='Binomial(n,'+str(nuclei)+','+str(prob_decay)+')')\n", "#ax1.scatter(unique,counts,marker='o',c='g',label='data')\n", "ax1.plot(unique,counts,marker='o',c='g',label='data')\n", "ax1.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Caso semplice con prob_decay = 0.001, 10**4 nuclei
\n", "Nota: ci possono essere problemi di overflow nel calcolare $\\lambda^n$ e i fattoriali
\n", "La distribuzione binomiale fallisce sui fattoriali" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nuclei = 10000\n", "ntry = 10000\n", "prob_decay = 0.001\n", "PseudoExp = np.random.uniform(size=(ntry,nuclei))/(1.-prob_decay)\n", "\n", "# Select values larger than 1 by using the floor method\n", "PseudoExp1 = np.floor(PseudoExp)\n", "print(PseudoExp1)\n", "# Count the decay taking the sum over lines\n", "decays = np.sum(PseudoExp1,axis=1).astype(int)\n", "\n", "ave = decays.mean()\n", "std_dev = decays.std()\n", "print(\"average counts:\",ave)\n", "print(\"standard deviation:\",std_dev)\n", "print(\"standard deviation squared:\",std_dev**2)\n", "unique, counts = np.unique(decays, return_counts=True)\n", "total = np.sum(counts)\n", "counts = counts/total" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(3*np.floor(ave))\n", "average_counts = nuclei*prob_decay\n", "yP = PoissonVal(t,average_counts)\n", "#yB = BinomialVal(t,nuclei,prob_decay)\n", "# check normalization \n", "print('Poisson Normalization:',np.sum(yP))\n", "\n", "fig2, ax2 = plt.subplots(figsize=(12, 5))\n", "ax2.set_title('Poisson distribution vs data for $\\lambda$ ='+str(average_counts))\n", "ax2.set_xlabel('n')\n", "ax2.set_ylabel('P(n)')\n", "ax2.grid(axis='x')\n", "#leg1 = '$\\lambda$ ='+str(average_counts)\n", "ax2.plot(t,yP,marker='o',c='b', linestyle = '--',label='$\\lambda$ ='+str(average_counts))\n", "#ax2.plot(t,yB,marker='o',c='k', linestyle = '--',label='Binomial(n,'+str(nuclei)+','+str(prob_decay)+')')\n", "ax2.plot(unique,counts,marker='o',c='g',label='data')\n", "ax2.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Poisson vs Gauss\n", "Confrontiamo\n", "$$ P_\\lambda(n) = \\frac{\\lambda^n}{n!} e^{-\\lambda}$$\n", "con\n", "$$ \\frac{1}{\\sqrt{2 \\pi \\lambda}} \\int^{n+1/2}_{n+1/2}dx \\, e^{-\\frac{1}{2}\\left(\\frac{x-\\lambda}{\\sqrt\\lambda}\\right)^2} = \\int^{n+1/2}_{n+1/2}dx\\, G(\\lambda,\\sqrt\\lambda,x)$$\n", "e con\n", "$$G(\\lambda,\\sqrt\\lambda,n)$$\n", "dove $G(\\mu,\\sigma,x)$ è la distribuzione normale.
\n", "Utilizzando la error function:\n", "$$\\operatorname{erf}(z) = \\frac{2}{\\sqrt{\\pi}}\\int^{z}_{0}e^{-x^2}dx$$\n", "si ha\n", "$$ \\frac{1}{\\sqrt{2 \\pi \\lambda}} \\int^{n+1/2}_{n+1/2}dx \\, e^{-\\frac{1}{2}\\left(\\frac{x-\\lambda}{\\sqrt\\lambda}\\right)^2} = \n", "\\frac{1}{2} \\left( \n", "\\operatorname{erf}\\left(\\frac{n+1/2-\\lambda}{\\sqrt{2\\lambda}}\\right) - \n", "\\operatorname{erf}\\left(\\frac{n-1/2-\\lambda}{\\sqrt{2\\lambda}}\\right)\\right)$$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(3,18)\n", "lam = 10." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "yGint = 0.5*(scipy.special.erf((t+0.5-lam)/np.sqrt(2*lam))-scipy.special.erf((t-0.5-lam)/np.sqrt(2*lam)))\n", "yG = np.exp(-0.5*(t-lam)**2/lam)/np.sqrt(2*np.pi*lam)\n", "yP = PoissonVal(t,lam)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig3, ax3 = plt.subplots(figsize=(12, 5))\n", "ax3.set_title('Poisson distribution vs Gauss distribution for $\\lambda$ = 10')\n", "ax3.set_xlabel('n')\n", "ax3.set_ylabel('P')\n", "ax3.grid(axis='x')\n", "ax3.plot(t,yP,marker='o',c='b', linestyle = '--',label='Poisson $\\lambda$ =10')\n", "ax3.plot(t,yGint,marker='o',c='r', linestyle = '--',label='Gauss_int $\\mu = 10$, $\\sigma = \\sqrt{10}$')\n", "ax3.plot(t,yG,marker='o',c='g', linestyle = '--',label='Gauss $\\mu = 10$, $\\sigma = \\sqrt{10}$')\n", "ax3.legend();" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(10,31)\n", "lam = 20.\n", "\n", "yGint = 0.5*(scipy.special.erf((t+0.5-lam)/np.sqrt(2*lam))-scipy.special.erf((t-0.5-lam)/np.sqrt(2*lam)))\n", "yP = PoissonVal(t,lam)\n", "\n", "fig3, ax3 = plt.subplots(figsize=(12, 5))\n", "ax3.set_title('Poisson distribution vs Gauss distribution for $\\lambda$ = 20')\n", "ax3.set_xlabel('n')\n", "ax3.set_ylabel('P')\n", "ax3.grid(axis='x')\n", "ax3.plot(t,yP,marker='o',c='b', linestyle = '--',label='Poisson $\\lambda$ =20')\n", "ax3.plot(t,yGint,marker='o',c='r', linestyle = '--',label='Gauss $\\mu = 20$, $\\sigma = \\sqrt{20}$')\n", "ax3.legend();" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(20,41)\n", "lam = 30.\n", "\n", "yGint = 0.5*(scipy.special.erf((t+0.5-lam)/np.sqrt(2*lam))-scipy.special.erf((t-0.5-lam)/np.sqrt(2*lam)))\n", "yP = PoissonVal(t,lam)\n", "\n", "fig3, ax3 = plt.subplots(figsize=(12, 5))\n", "ax3.set_title('Poisson distribution vs Gauss distribution for $\\lambda$ = 30')\n", "ax3.set_xlabel('n')\n", "ax3.set_ylabel('P')\n", "ax3.grid(axis='x')\n", "ax3.plot(t,yP,marker='o',c='b', linestyle = '--',label='Poisson $\\lambda$ =30')\n", "ax3.plot(t,yGint,marker='o',c='r', linestyle = '--',label='Gauss $\\mu = 30$, $\\sigma = \\sqrt{30}$')\n", "ax3.legend();" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(80,121)\n", "lam = 100.\n", "\n", "yGint = 0.5*(scipy.special.erf((t+0.5-lam)/np.sqrt(2*lam))-scipy.special.erf((t-0.5-lam)/np.sqrt(2*lam)))\n", "yP = PoissonVal(t,lam)\n", "\n", "fig3, ax3 = plt.subplots(figsize=(12, 5))\n", "ax3.set_title('Poisson distribution vs Gauss distribution for $\\lambda$ = 100')\n", "ax3.set_xlabel('n')\n", "ax3.set_ylabel('P')\n", "ax3.grid(axis='x')\n", "ax3.plot(t,yP,marker='o',c='b', linestyle = '--',label='Poisson $\\lambda$ =100')\n", "ax3.plot(t,yGint,marker='o',c='r', linestyle = '--',label='Gauss $\\mu = 100$, $\\sigma = \\sqrt{100}$')\n", "ax3.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test computing factorials\n", "Calcolare i termini della distribuzione binomiale \n", "$$B_p(n,N) = \\frac{N!}{n!(N-n)!}\\, p^n (1-p)^{N-n}$$ \n", "per grandi n pone dei problemi. Anche se gli interi in Python sono arbitrariamente grandi, devono essere convertiti in float per completare il calcolo." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bigFact = scipy.special.factorial(1000)\n", "bigFact" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### float e longdouble non sono sufficienti a contenere bigFact" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.float(bigFact)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.longdouble(bigFact)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scipy.special.factorial(10)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.asarray(scipy.special.factorial(10),dtype=np.int64)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.asarray(scipy.special.factorial(1000),dtype=np.int64)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def my_fact(n):\n", " res=1\n", " for i in range(2,n+1):\n", " res *=i\n", " return res\n", "\n", "my_fact(10)\n", "# my_fact(10000) # funziona" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Passiamo ai logaritmi" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import math\n", "n1 = my_fact(10000);\n", "n2 = my_fact(5000);\n", "res = math.log(n1)-2.*math.log(n2)\n", "res" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res1 = res +2*5000*math.log(0.5)\n", "res1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res1 = math.exp(res1)\n", "res1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### test con n1 = 200 (con n1 = 300 non riesco più a calcolare senza logaritmi)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = math.log(n1f)-2.*math.log(n2f)+n2*math.log(p)+n2*math.log(1.-p)\n", "res = math.exp(res)\n", "print(res)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Torniamo a Poisson vs Binomial vs data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def BinomialVal2(n,ntot,p):\n", " '''\n", " Returns value of the Binomial distribution for n successes in ntot tries with\n", " probability p for a success.\n", " '''\n", " n1 = ntot-n\n", " logres = math.log(my_fact(ntot)) - math.log(my_fact(n)) \\\n", " - math.log(my_fact(n1)) + n*math.log(p) + n1*math.log(1.-p)\n", " return math.exp(logres) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nuclei = 10000\n", "ntry = 10000\n", "prob_decay = 0.001\n", "PseudoExp = np.random.uniform(size=(ntry,nuclei))/(1.-prob_decay)\n", "\n", "# Select values larger than 1 by using the floor method\n", "PseudoExp1 = np.floor(PseudoExp)\n", "print(PseudoExp1)\n", "# Count the decay taking the sum over lines\n", "decays = np.sum(PseudoExp1,axis=1).astype(int)\n", "\n", "ave = decays.mean()\n", "std_dev = decays.std()\n", "print(\"average counts:\",ave)\n", "print(\"standard deviation:\",std_dev)\n", "print(\"standard deviation squared:\",std_dev**2)\n", "unique, counts = np.unique(decays, return_counts=True)\n", "total = np.sum(counts)\n", "counts = counts/total" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.arange(3*np.floor(ave))\n", "average_counts = nuclei*prob_decay\n", "yP = PoissonVal(t,average_counts)\n", "t1 = range(3*math.floor(ave))\n", "yB = list([BinomialVal2(x,nuclei,prob_decay) for x in t1])\n", "\n", "# check normalization \n", "print('Poisson Normalization:',np.sum(yP))\n", "\n", "fig2, ax2 = plt.subplots(figsize=(12, 5))\n", "ax2.set_title('Poisson distribution vs data for $\\lambda$ ='+str(average_counts))\n", "ax2.set_xlabel('n')\n", "ax2.set_ylabel('P(n)')\n", "ax2.grid(axis='x')\n", "#leg1 = '$\\lambda$ ='+str(average_counts)\n", "ax2.plot(t,yP,marker='o',c='b', linestyle = '--',label='$\\lambda$ ='+str(average_counts))\n", "ax2.plot(t1,yB,marker='x',markersize=15, c='k', linestyle = '--',label='Binomial(n,'+str(nuclei)+','+str(prob_decay)+')')\n", "ax2.plot(unique,counts,marker='o', c='g',label='data')\n", "ax2.legend();" ] }, { "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.9.12" } }, "nbformat": 4, "nbformat_minor": 4 }