From 8425ff90daae87356e88daa26c9a316c9b9309eb Mon Sep 17 00:00:00 2001 From: Jacopo Basso Ricci Date: Fri, 6 Nov 2026 10:26:45 +0100 Subject: [PATCH] Aggiunta dati e versione finale di focale lente --- csv/focaleLente/lenti per csv 2(252B_m1).csv | 11 + csv/focaleLente/lenti per csv 2(252B_m2).csv | 11 + csv/focaleLente/lenti per csv 2(252B_m3).csv | 11 + csv/focaleLente/lenti per csv 2(48B_m1).csv | 8 + csv/focaleLente/lenti per csv 2(48B_m2).csv | 9 + csv/focaleLente/lenti per csv 2(48B_m3).csv | 11 + .../lenti per csv 2(dimensione freccia).csv | 2 + .../lenti per csv 2(spessore 252B).csv | 2 + .../lenti per csv 2(spessore 48B).csv | 2 + focaleLente.py | 849 +++++++++++------- pyproject.toml | 1 + statistica.py | 126 ++- 12 files changed, 715 insertions(+), 328 deletions(-) create mode 100644 csv/focaleLente/lenti per csv 2(252B_m1).csv create mode 100644 csv/focaleLente/lenti per csv 2(252B_m2).csv create mode 100644 csv/focaleLente/lenti per csv 2(252B_m3).csv create mode 100644 csv/focaleLente/lenti per csv 2(48B_m1).csv create mode 100644 csv/focaleLente/lenti per csv 2(48B_m2).csv create mode 100644 csv/focaleLente/lenti per csv 2(48B_m3).csv create mode 100644 csv/focaleLente/lenti per csv 2(dimensione freccia).csv create mode 100644 csv/focaleLente/lenti per csv 2(spessore 252B).csv create mode 100644 csv/focaleLente/lenti per csv 2(spessore 48B).csv diff --git a/csv/focaleLente/lenti per csv 2(252B_m1).csv b/csv/focaleLente/lenti per csv 2(252B_m1).csv new file mode 100644 index 0000000..1baa368 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(252B_m1).csv @@ -0,0 +1,11 @@ +p ;q +38,8;68,6 +66,3;41,0 +42,6;63,6 +49,4;52,9 +64,7;42,2 +48,7;54,1 +51,1;50,1 +48,7;55,0 +50,7;51,0 +63,9;42,2 diff --git a/csv/focaleLente/lenti per csv 2(252B_m2).csv b/csv/focaleLente/lenti per csv 2(252B_m2).csv new file mode 100644 index 0000000..3fcdb57 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(252B_m2).csv @@ -0,0 +1,11 @@ +q ;h +68,6;21,70 +41,0;8,20 +63,6;18,90 +41,4;8,00 +42,7;8,75 +51,8;13,00 +55,4;15,05 +43,2;9,00 +39,6;6,60 +41,2;7,30 diff --git a/csv/focaleLente/lenti per csv 2(252B_m3).csv b/csv/focaleLente/lenti per csv 2(252B_m3).csv new file mode 100644 index 0000000..2fed514 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(252B_m3).csv @@ -0,0 +1,11 @@ +a ;d1 ;d2 +108,7;40,0;68,0 +106,6;41,0;64,0 +103,7;42,0;57,5 +110,3;44,5;55,5 +103,2;43,5;61,0 +107,0;42,5;67,0 +110,5;38,5;70,0 +104,0;45,0;59,0 +101,7;45,0;56,0 +105,0;42,0;58,0 diff --git a/csv/focaleLente/lenti per csv 2(48B_m1).csv b/csv/focaleLente/lenti per csv 2(48B_m1).csv new file mode 100644 index 0000000..96bd898 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(48B_m1).csv @@ -0,0 +1,8 @@ +p;q +5,4;73,7 +5,7;50,5 +5,4;86,7 +5,2;98,8 +5,3;105,3 +5,4;77,1 +5,5;66,9 diff --git a/csv/focaleLente/lenti per csv 2(48B_m2).csv b/csv/focaleLente/lenti per csv 2(48B_m2).csv new file mode 100644 index 0000000..fa9f780 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(48B_m2).csv @@ -0,0 +1,9 @@ +q;h +11,1;22,00 +20,3;42,20 +27,9;62,45 +36,2;85,10 +43,9;106,15 +43,7;105,00 +48,5;105,50 +7,3;8,30 diff --git a/csv/focaleLente/lenti per csv 2(48B_m3).csv b/csv/focaleLente/lenti per csv 2(48B_m3).csv new file mode 100644 index 0000000..513cf12 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(48B_m3).csv @@ -0,0 +1,11 @@ +a ;d1 ;d2 +24,6;18,4;6,7 +26,6;20,8;6,3 +32,2;26,5;6,0 +37,9;32,6;5,4 +45,4;40,5;5,3 +55,2;49,8;5,2 +39,5;34,1;5,4 +33,2;27,2;5,6 +23,7;17,0;6,6 +41,2;35,7;5,5 \ No newline at end of file diff --git a/csv/focaleLente/lenti per csv 2(dimensione freccia).csv b/csv/focaleLente/lenti per csv 2(dimensione freccia).csv new file mode 100644 index 0000000..5bc8c41 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(dimensione freccia).csv @@ -0,0 +1,2 @@ +freccia1;freccia2;freccia3;freccia4;freccia5;freccia6;freccia7;freccia8;freccia9;freccia10;freccia11 +11,90;11,70;11,80;11,95;11,85;11,95;11,90;11,85;11,90;11,95;11,85 diff --git a/csv/focaleLente/lenti per csv 2(spessore 252B).csv b/csv/focaleLente/lenti per csv 2(spessore 252B).csv new file mode 100644 index 0000000..bf89084 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(spessore 252B).csv @@ -0,0 +1,2 @@ +s252B1;s252B2;s252B3;s252B4;s252B5;s252B6;s252B7;s252B8 +6,45;6,35;6,40;6,45;6,50;6,45;6,40;6,45 diff --git a/csv/focaleLente/lenti per csv 2(spessore 48B).csv b/csv/focaleLente/lenti per csv 2(spessore 48B).csv new file mode 100644 index 0000000..9445ed1 --- /dev/null +++ b/csv/focaleLente/lenti per csv 2(spessore 48B).csv @@ -0,0 +1,2 @@ +s48B1;s48B2;s48B3;s48B4;s48B5;s48B6;s48B7;s48B8 +6,45;6,45;6,50;6,45;6,40;6,40;6,45;6,40 diff --git a/focaleLente.py b/focaleLente.py index 77f0167..d0c53ce 100644 --- a/focaleLente.py +++ b/focaleLente.py @@ -1,402 +1,627 @@ import marimo -__generated_with = "0.25.0" -app = marimo.App(width="medium") +__generated_with = "0.25.2" +app = marimo.App(width="medium", auto_download=["html"]) @app.cell def _(): import marimo as mo - - return (mo,) - - -@app.cell(hide_code=True) -def _(mo): - mo.md(r""" - # Misura della distanza focale di una lente convergente - - Quattro stime, tutte basate sull'equazione della lente sottile $\frac{1}{f}=\frac{1}{p}+\frac{1}{q}$ e sull'ingrandimento $m=\frac{q}{p}=\frac{l_{Freccia}}{L}$ (con $L$ = lunghezza dell'oggetto): - - 1. **Punti coniugati**: $f=\dfrac{p\,q}{p+q}$ - 2. **a) $p$ e ingrandimento**: $f=\dfrac{p\,m}{1+m}$ - 2. **b) $q$ e ingrandimento**: $f=\dfrac{q}{1+m}$ - 3. **Bessel**: $f=\dfrac{D^2-d^2}{4D}$ (vedi sotto) - - Le incertezze sono propagate con `Data.calc_var` di `statistica.py`. - """) - return - - -@app.cell -def _(): import numpy as np import pandas as pd import sympy as sp from scipy import stats + from pathlib import Path + import statistica - # Esegue statistica.py _outputs, defs = statistica.app.run() - Data = defs["Data"] compat = defs["compat"] + media_pesata = defs["media_pesata"] + media_pesata_dersimonian_laird = defs["media_pesata_dersimonian_laird"] + rimuovi_outlier = defs["rimuovi_outlier"] plot_gauss = defs["plot_gauss"] - return Data, compat, np, pd, plot_gauss, sp, stats + cartella_csv = Path(__file__).resolve().parent / "csv" / "focaleLente" + cartella_immagini = Path(__file__).resolve().parent / "LaTeX" / "immagini" + return ( + Data, + cartella_csv, + cartella_immagini, + compat, + media_pesata, + mo, + np, + pd, + plot_gauss, + rimuovi_outlier, + sp, + ) @app.cell(hide_code=True) def _(mo): mo.md(r""" - ## Dati + # Misura della distanza focale - **Le righe sono a coppie (0-1, 2-3, 4-5, ...)** per il metodo di Bessel: le due righe di una coppia hanno la **stessa distanza sorgente-schermo** $D=p+q$ e corrispondono alle due posizioni della lente con immagine nitida. + Sono usati solo i CSV presenti in `csv/focaleLente`. - In una posizione la lente è vicina alla sorgente (immagine ingrandita, $p$ piccolo), nell'altra è vicina allo schermo (immagine rimpicciolita, $p$ grande): le due posizioni sono **coniugate**, cioè $p$ e $q$ si scambiano ($p_B=q_A$, $q_B=p_A$) e gli ingrandimenti sono uno il reciproco dell'altro. + - **m1, punti coniugati:** $f = pq/(p+q)$ + - **m2, $q$ e ingrandimento:** $f = q/(1+h/L)$ + - **m3, Bessel:** $f = d_1d_2/a$ - Ordine delle righe dentro la coppia indifferente. + Il metodo con $p$ e ingrandimento non viene usato perché non è stata + misurata la grandezza necessaria. + + Le distanze da metro sono in cm con errore strumentale `0.1 cm`. + Le misure da calibro ventesimale sono in mm con errore strumentale + `0.05 mm`. """) return @app.cell -def _(Data, pd): - dfFenditura = pd.DataFrame({ - "l1": [10.05,], - "l2": [10.00,], - "l3": [10.03,], - "l4": [9.94, ], - "l5": [9.99, ], - }) +def _(cartella_csv, pd): + def leggi_csv(nome): + dati = pd.read_csv(cartella_csv / nome, sep=";", decimal=",") + dati.columns = dati.columns.astype(str).str.strip() + dati = dati.dropna(axis=1, how="all") + return dati.apply(pd.to_numeric, errors="coerce") - datiFenditura = Data(dfFenditura) - datiFenditura.stat_chauv(["l"], [0.02]) - datiFenditura - - lunghezzaOggetto = datiFenditura.analisi_stat.l[0] # variabile che contiene la lunghezza dell'oggetto - ulunghezzaOggetto = datiFenditura.analisi_stat.ul[0] # variabile che contiene l'errore sulla lunghezza dell'oggetto - return lunghezzaOggetto, ulunghezzaOggetto + m1_48B = leggi_csv("lenti per csv 2(48B_m1).csv") + m1_252B = leggi_csv("lenti per csv 2(252B_m1).csv") + m2_48B = leggi_csv("lenti per csv 2(48B_m2).csv") + m2_252B = leggi_csv("lenti per csv 2(252B_m2).csv") + m3_48B = leggi_csv("lenti per csv 2(48B_m3).csv") + m3_252B = leggi_csv("lenti per csv 2(252B_m3).csv") + spessore_48B_df = leggi_csv("lenti per csv 2(spessore 48B).csv") + spessore_252B_df = leggi_csv("lenti per csv 2(spessore 252B).csv") + dimensione_freccia_df = leggi_csv("lenti per csv 2(dimensione freccia).csv") + return ( + dimensione_freccia_df, + m1_252B, + m1_48B, + m2_252B, + m2_48B, + m3_252B, + m3_48B, + spessore_252B_df, + spessore_48B_df, + ) @app.cell -def _(pd): - df = pd.DataFrame({ - "p1": [23.82, 41.22, 20.58, 53.91, 19.38, 65.51], - "p2": [23.39, 41.62, 20.94, 54.35, 19.07, 65.62], - "p3": [23.49, 41.54, 20.40, 54.72, 19.43, 65.70], - "p4": [23.57, 41.57, 20.84, 54.32, 19.79, 65.51], - "p5": [23.33, 41.21, 20.32, 54.17, 19.37, 65.19], +def _(Data, dimensione_freccia_df, spessore_252B_df, spessore_48B_df): + _df_spessore_48B = Data(spessore_48B_df) + _df_spessore_48B.stat_chauv("s48B", 0.05, False) + spessore_48B = _df_spessore_48B.analisi_stat.s48B[0] /10 + uspessore_48B = _df_spessore_48B.analisi_stat.us48B[0]/10 - "q1": [41.51, 23.98, 54.07, 21.30, 65.28, 19.65], - "q2": [41.51, 23.53, 53.91, 20.80, 65.19, 19.72], - "q3": [40.99, 23.37, 54.71, 20.76, 65.23, 19.58], - "q4": [41.82, 24.09, 54.80, 20.81, 65.46, 19.17], - "q5": [41.69, 23.47, 54.17, 20.69, 65.19, 19.39], + _df_spessore_252B = Data(spessore_252B_df) + _df_spessore_252B.stat_chauv("s252B", 0.05, False) + spessore_252B = _df_spessore_252B.analisi_stat.s252B[0] /10 + uspessore_252B = _df_spessore_252B.analisi_stat.us252B[0]/10 - "lFreccia1": [17.64, 5.58, 26.22, 3.80, 33.77, 2.94], - "lFreccia2": [17.67, 5.64, 26.17, 3.75, 33.69, 2.95], - "lFreccia3": [17.70, 5.54, 26.21, 3.84, 33.70, 3.03], - "lFreccia4": [17.66, 5.71, 26.14, 3.81, 33.77, 3.04], - "lFreccia5": [17.66, 5.64, 26.09, 3.88, 33.77, 3.00], - }) + _df_freccia = Data(dimensione_freccia_df) + _df_freccia.stat_chauv("freccia", 0.05, False) + freccia = _df_freccia.analisi_stat.freccia[0] /10 + ufreccia = _df_freccia.analisi_stat.ufreccia[0]/10 + return ( + freccia, + spessore_252B, + spessore_48B, + ufreccia, + uspessore_252B, + uspessore_48B, + ) - # incertezze strumentali (risoluzione) di p, q, lFreccia - err_strum = {"p": 0.1, "q": 0.1, "lFreccia": 0.05} - df - return df, err_strum + +@app.cell(hide_code=True) +def _( + freccia, + mo, + spessore_252B, + spessore_48B, + ufreccia, + uspessore_252B, + uspessore_48B, +): + mo.md("## Misure accessorie") + (spessore_48B, uspessore_48B), (spessore_252B, uspessore_252B), (freccia, ufreccia) + return + + +@app.cell +def _( + Data, + m1_252B, + m1_48B, + spessore_252B, + spessore_48B, + uspessore_252B, + uspessore_48B, +): + data48Bm1 = Data(m1_48B.copy()) + data48Bm1.analisi_stat = m1_48B.copy() + data48Bm1.analisi_stat["up"] = 0.1 + data48Bm1.analisi_stat["uq"] = 1 #Potevamo muovere lo schermo di ... e non notavamo differenze + data48Bm1.analisi_stat["s"] = spessore_48B + data48Bm1.analisi_stat["us"] = uspessore_48B + + data252Bm1 = Data(m1_252B.copy()) + data252Bm1.analisi_stat = m1_252B.copy() + data252Bm1.analisi_stat["up"] = 1 #Potevamo muovere lo schermo di ... e non notavamo differenze + data252Bm1.analisi_stat["uq"] = 1 #Potevamo muovere lo schermo di ... e non notavamo differenze + data252Bm1.analisi_stat["s"] = spessore_252B + data252Bm1.analisi_stat["us"] = uspessore_252B + + data48Bm1.analisi_stat, data252Bm1.analisi_stat + return data252Bm1, data48Bm1 + + +@app.cell +def _(sp): + p, q, s = sp.symbols("p q s", positive=True) + peff, qeff = sp.symbols("peff qeff", positive=True) + + peff = p - (s/2) + qeff = q + (s/2) + return peff, qeff + + +@app.cell +def _(data252Bm1, data48Bm1, peff, qeff): + f1 = peff * qeff / (peff + qeff) + + data48Bm1.calc_var(peff, "peff") + data252Bm1.calc_var(peff, "peff") + data48Bm1.calc_var(qeff, "qeff") + data252Bm1.calc_var(qeff, "qeff") + data48Bm1.calc_var(f1, "f1") + data252Bm1.calc_var(f1, "f1") + + risultato48Bm1 = data48Bm1.analisi_stat[["peff", "upeff", "qeff", "uqeff", "f1", "uf1"]] + risultato252Bm1 = data252Bm1.analisi_stat[["peff", "upeff", "qeff", "uqeff", "f1", "uf1"]] + return risultato252Bm1, risultato48Bm1 @app.cell(hide_code=True) def _(mo): mo.md(r""" - ## Analisi statistica (criterio di Chauvenet) + ## Metodo 1: punti coniugati - Per ogni configurazione: media, incertezza statistica + strumentale, ed - eventuali outlier scartati. Alla tabella si aggiunge la lunghezza - dell'oggetto, comune a tutte le righe. + `p` e `q` sono in cm, quindi anche `f1` è in cm. """) return @app.cell -def _(Data, df, err_strum, lunghezzaOggetto, ulunghezzaOggetto): - dati = Data(df) - dati.stat_chauv(list(err_strum.keys()), list(err_strum.values())) +def _(risultato252Bm1, risultato48Bm1): + risultato48Bm1, risultato252Bm1 + return - dati.analisi_stat["lunghezzaOggetto"] = lunghezzaOggetto - dati.analisi_stat["ulunghezzaOggetto"] = ulunghezzaOggetto - dati.analisi_stat - return (dati,) +@app.cell +def _( + Data, + freccia, + m2_252B, + m2_48B, + spessore_252B, + spessore_48B, + ufreccia, + uspessore_252B, + uspessore_48B, +): + data48Bm2 = Data(m2_48B.copy()) + data48Bm2.analisi_stat = m2_48B.copy() + data48Bm2.analisi_stat["uq"] = 1 #La freccia era sempre a fuoco per ogni q + data48Bm2.analisi_stat["h"] /= 10 + data48Bm2.analisi_stat["uh"] = 0.1 #No guide fisse e calibro nell'vuoto siderale + data48Bm2.analisi_stat["L"] = freccia + data48Bm2.analisi_stat["uL"] = ufreccia + data48Bm2.analisi_stat["s"] = spessore_48B + data48Bm2.analisi_stat["us"] = uspessore_48B + + data252Bm2 = Data(m2_252B.copy()) + data252Bm2.analisi_stat = m2_252B.copy() + data252Bm2.analisi_stat["uq"] = 0.1 + data252Bm2.analisi_stat["h"] /= 10 + data252Bm2.analisi_stat["uh"] = 0.02 + data252Bm2.analisi_stat["L"] = freccia + data252Bm2.analisi_stat["uL"] = ufreccia + data252Bm2.analisi_stat["s"] = spessore_252B + data252Bm2.analisi_stat["us"] = uspessore_252B + return data252Bm2, data48Bm2 + + +@app.cell +def _(data252Bm2, data48Bm2, np, qeff, sp): + h, L, theta= sp.symbols("h L theta", positive=True) + + data48Bm2.calc_var(qeff, "qeff") + data252Bm2.calc_var(qeff, "qeff") + + data48Bm2.analisi_stat["theta"] = np.radians(15) + data48Bm2.analisi_stat["utheta"] = (np.radians(15) - np.radians(5)) / 2 # data48Bm2.analisi_stat["theta"] / (2* np.sqrt(3)) + + data252Bm2.analisi_stat["theta"] = np.radians(15) + data252Bm2.analisi_stat["utheta"] = (np.radians(15) - np.radians(5)) / 2 + + + f2bCorretto = qeff / (1 + h * sp.cos(theta) / L) # / 1.0353 # [1; 1.0353] + f2b = qeff / (1 + h / L) #Senza correzione si theta per vedere come è + + data48Bm2.calc_var(f2b, "f2b_senza") + data252Bm2.calc_var(f2b, "f2b_senza") + data48Bm2.calc_var(f2bCorretto, "f2b") + data252Bm2.calc_var(f2bCorretto, "f2b") + + risultato48Bm2 = data48Bm2.analisi_stat[["qeff", "uqeff", "h", "uh", "f2b", "uf2b"]] + risultato252Bm2 = data252Bm2.analisi_stat[["qeff", "uqeff", "h", "uh", "f2b", "uf2b"]] + risultato48Bm2_senza = data48Bm2.analisi_stat[["qeff", "uqeff", "h", "uh", "f2b_senza", "uf2b_senza"]] + risultato252Bm2_senza = data252Bm2.analisi_stat[["qeff", "uqeff", "h", "uh", "f2b_senza", "uf2b_senza"]] + + data48Bm2.analisi_stat, data252Bm2.analisi_stat, + return ( + risultato252Bm2, + risultato252Bm2_senza, + risultato48Bm2, + risultato48Bm2_senza, + ) @app.cell(hide_code=True) def _(mo): mo.md(r""" - ## Metodo 1 — punti coniugati - """) - return + ## Metodo 2: $q$ e ingrandimento - -@app.cell(hide_code=True) -def _(f1, mo, sp): - mo.md(f""" - $f = {sp.latex(f1)}$ + `q` e `L` sono in cm e `h` è in mm. Il rapporto `h/L` è adimensionale, + perché confronta la dimensione dell'immagine con quella dell'oggetto. + Il risultato `f2b` è in cm. """) return @app.cell -def _(dati, sp): - _p, _q = sp.symbols("p q", positive=True) +def _(risultato252Bm2, risultato48Bm2): + risultato48Bm2, risultato252Bm2 + return - f1 = _p * _q / (_p + _q) - dati.calc_var(f1, "f1") +@app.cell +def _(Data, m3_252B, m3_48B): + data48Bm3 = Data(m3_48B.copy()) + data48Bm3.analisi_stat = m3_48B.copy() + data48Bm3.analisi_stat["ua"] = 0.3 + data48Bm3.analisi_stat["ud1"] = 0.1 # Non notavamo differenze in 0.3 cm + data48Bm3.analisi_stat["ud2"] = 0.1 - res_1 = dati.analisi_stat[["f1", "uf1", "%uf1_p", "%uf1_q"]] - res_1 - return f1, res_1 + data252Bm3 = Data(m3_252B.copy()) + data252Bm3.analisi_stat = m3_252B.copy() + data252Bm3.analisi_stat["ua"] = 0.3 + data252Bm3.analisi_stat["ud1"] = 1 # No differenze notate in 1 cm + data252Bm3.analisi_stat["ud2"] = 1 + return data252Bm3, data48Bm3 + + +@app.cell +def _(data252Bm3, data48Bm3, sp): + a, d1, d2 = sp.symbols("a d1 d2", positive=True, real=True) + delta = d1 - d2 + + data48Bm3.calc_var(delta, "delta") + data252Bm3.calc_var(delta, "delta") + + f3 = (a**2 - delta**2) / (4*a) + + data48Bm3.calc_var(f3, "f3") + data252Bm3.calc_var(f3, "f3") + + risultato48Bm3 = data48Bm3.analisi_stat[["a", "ua", "d1", "ud1", "d2", "ud2", "f3", "uf3"]] + risultato252Bm3 = data252Bm3.analisi_stat[["a", "ua", "d1", "ud1", "d2", "ud2", "f3", "uf3"]] + + data48Bm3.analisi_stat + return risultato252Bm3, risultato48Bm3 @app.cell(hide_code=True) def _(mo): mo.md(r""" - ## Metodo 2a — $p$ e ingrandimento + ## Metodo 3: Bessel + + `a`, `d1` e `d2` sono in cm. Si usa direttamente """) return +@app.cell +def _(risultato252Bm3, risultato48Bm3): + risultato48Bm3, risultato252Bm3 + return + + +@app.cell +def _( + rimuovi_outlier, + risultato252Bm1, + risultato252Bm2, + risultato252Bm3, + risultato48Bm1, + risultato48Bm2, + risultato48Bm3, +): + f48B_m1_puliti, out48B_m1 = rimuovi_outlier( + risultato48Bm1["f1"].to_numpy(), risultato48Bm1["uf1"].mean() + ) + f48B_m2_puliti, out48B_m2 = rimuovi_outlier( + risultato48Bm2["f2b"].to_numpy(), risultato48Bm2["uf2b"].mean() + ) + f48B_m3_puliti, out48B_m3 = rimuovi_outlier( + risultato48Bm3["f3"].to_numpy(), risultato48Bm3["uf3"].mean() + ) + + f252B_m1_puliti, out252B_m1 = rimuovi_outlier( + risultato252Bm1["f1"].to_numpy(), risultato252Bm1["uf1"].mean() + ) + f252B_m2_puliti, out252B_m2 = rimuovi_outlier( + risultato252Bm2["f2b"].to_numpy(), risultato252Bm2["uf2b"].mean() + ) + f252B_m3_puliti, out252B_m3 = rimuovi_outlier( + risultato252Bm3["f3"].to_numpy(), risultato252Bm3["uf3"].mean() + ) + + outlier_fuochi = { + "48B m1": out48B_m1, + "48B m2": out48B_m2, + "48B m3": out48B_m3, + "252B m1": out252B_m1, + "252B m2": out252B_m2, + "252B m3": out252B_m3, + } + return ( + f252B_m1_puliti, + f252B_m2_puliti, + f252B_m3_puliti, + f48B_m1_puliti, + f48B_m2_puliti, + f48B_m3_puliti, + outlier_fuochi, + ) + + @app.cell(hide_code=True) -def _(f2a, mo, sp): - mo.md(f""" - $f = {sp.latex(f2a)}$ - """) +def _(mo, outlier_fuochi): + mo.md("## Outlier dei fuochi") + outlier_fuochi return @app.cell -def _(dati, sp): - _p, _lF, _L = sp.symbols("p lFreccia lunghezzaOggetto", positive=True) - _m = _lF / _L +def _( + f252B_m1_puliti, + f252B_m2_puliti, + f252B_m3_puliti, + f48B_m1_puliti, + f48B_m2_puliti, + f48B_m3_puliti, + media_pesata, + np, + pd, + risultato252Bm1, + risultato252Bm2, + risultato252Bm3, + risultato48Bm1, + risultato48Bm2, + risultato48Bm3, +): + f48B_m1, u48B_m1, chi2_48B_m1, dof_48B_m1, P_48B_m1 = media_pesata( + f48B_m1_puliti, risultato48Bm1["uf1"].loc[risultato48Bm1["f1"].isin(f48B_m1_puliti)] + ) + f48B_m2, u48B_m2, chi2_48B_m2, dof_48B_m2, P_48B_m2 = media_pesata( + f48B_m2_puliti, risultato48Bm2["uf2b"].loc[risultato48Bm2["f2b"].isin(f48B_m2_puliti)] + ) + f48B_m3, u48B_m3, chi2_48B_m3, dof_48B_m3, P_48B_m3 = media_pesata( + f48B_m3_puliti, risultato48Bm3["uf3"].loc[risultato48Bm3["f3"].isin(f48B_m3_puliti)] + ) - f2a = _p * _m / (1 + _m) + f252B_m1, u252B_m1, chi2_252B_m1, dof_252B_m1, P_252B_m1 = media_pesata( + f252B_m1_puliti, risultato252Bm1["uf1"].loc[risultato252Bm1["f1"].isin(f252B_m1_puliti)] + ) + f252B_m2, u252B_m2, chi2_252B_m2, dof_252B_m2, P_252B_m2 = media_pesata( + f252B_m2_puliti, risultato252Bm2["uf2b"].loc[risultato252Bm2["f2b"].isin(f252B_m2_puliti)] + ) + f252B_m3, u252B_m3, chi2_252B_m3, dof_252B_m3, P_252B_m3 = media_pesata( + f252B_m3_puliti, risultato252Bm3["uf3"].loc[risultato252Bm3["f3"].isin(f252B_m3_puliti)] + ) - dati.calc_var(f2a, "f2a") + risultati_finali = pd.DataFrame({ + "48B - m1": [f48B_m1, u48B_m1, chi2_48B_m1, dof_48B_m1, P_48B_m1], + "48B - m2": [f48B_m2, u48B_m2, chi2_48B_m2, dof_48B_m2, P_48B_m2], + "48B - m3": [f48B_m3, u48B_m3, chi2_48B_m3, dof_48B_m3, P_48B_m3], + "252B - m1": [f252B_m1, u252B_m1, chi2_252B_m1, dof_252B_m1, P_252B_m1], + "252B - m2": [f252B_m2, u252B_m2, chi2_252B_m2, dof_252B_m2, P_252B_m2], + "252B - m3": [f252B_m3, u252B_m3, chi2_252B_m3, dof_252B_m3, P_252B_m3], + }, index=["f [cm]", "u [cm]", "chi2", "dof", "P"]).T - res_2a = dati.analisi_stat[ - ["f2a", "uf2a", "%uf2a_p", "%uf2a_lFreccia", "%uf2a_lunghezzaOggetto"] - ] - res_2a - return f2a, res_2a + risultati_finali["u_corr [cm]"] = risultati_finali["u [cm]"] * np.sqrt( + np.maximum(1, risultati_finali["chi2"] / risultati_finali["dof"]) + ) + risultati_finali + return (risultati_finali,) @app.cell(hide_code=True) -def _(mo): - mo.md(r""" - ## Metodo 2b — $q$ e ingrandimento - """) - return +def _( + cartella_immagini, + media_pesata, + plot_gauss, + rimuovi_outlier, + risultati_finali, + risultato252Bm2_senza, + risultato48Bm2_senza, +): + f48B_m2_senza_puliti, _ = rimuovi_outlier( + risultato48Bm2_senza["f2b_senza"].to_numpy(), + risultato48Bm2_senza["uf2b_senza"].mean(), + ) + f252B_m2_senza_puliti, _ = rimuovi_outlier( + risultato252Bm2_senza["f2b_senza"].to_numpy(), + risultato252Bm2_senza["uf2b_senza"].mean(), + ) + f48B_m2_senza, u48B_m2_senza, _, _, _ = media_pesata( + f48B_m2_senza_puliti, + risultato48Bm2_senza["uf2b_senza"].loc[ + risultato48Bm2_senza["f2b_senza"].isin(f48B_m2_senza_puliti) + ], + ) + f252B_m2_senza, u252B_m2_senza, _, _, _ = media_pesata( + f252B_m2_senza_puliti, + risultato252Bm2_senza["uf2b_senza"].loc[ + risultato252Bm2_senza["f2b_senza"].isin(f252B_m2_senza_puliti) + ], + ) - -@app.cell(hide_code=True) -def _(f2b, mo, sp): - mo.md(f""" - $f = {sp.latex(f2b)}$ - """) - return - - -@app.cell -def _(dati, sp): - _q, _lF, _L = sp.symbols("q lFreccia lunghezzaOggetto", positive=True) - _m = _lF / _L - - f2b = _q / (1 + _m) - - dati.calc_var(f2b, "f2b") - - res_2b = dati.analisi_stat[ - ["f2b", "uf2b", "%uf2b_q", "%uf2b_lFreccia", "%uf2b_lunghezzaOggetto"] - ] - res_2b - return f2b, res_2b - - -@app.cell(hide_code=True) -def _(mo): - mo.md(r""" - ## Metodo 3 — Bessel - - Sorgente e schermo sono a distanza fissa $D>4f$. Esistono due posizioni della lente, distanti $d$ tra loro, con immagine nitida. - - Per la prima (A) $p_A=\frac{D-d}{2}$, $q_A=\frac{D+d}{2}$, quindi - $$f=\frac{p_Aq_A}{p_A+q_A}=\frac{D^2-d^2}{4D}$$ - - Ogni coppia di righe dà **una** stima di $f$. - Le quattro misure $p_A,q_A,p_B,q_B$ misurano tutte $D$ e $d$ (con A = posizione a $p$ minore): - - $$D=\frac{p_A+q_A+p_B+q_B}{2},\qquad d=\frac{p_B-p_A+q_A-q_B}{2}$$ - - **Attenzione alla propagazione:** $D$ e $d$ sono combinazioni delle stesse quattro misure, quindi *non* sono indipendenti. - - Propagare $f(D,d)$ trattandoli come indipendenti sbaglierebbe l'incertezza: $f$ viene perciò propagata direttamente da $p_A,q_A,p_B,q_B$ (che sono misure indipendenti). $D$ e $d$ (con le loro incertezze) sono calcolati solo per i controlli. - """) - return - - -@app.cell -def _(dati, pd): - _s = dati.analisi_stat - assert len(_s) % 2 == 0, "Bessel: servono righe a coppie (0-1, 2-3, ...)" - - _righe, _etichette = [], [] - for _k in range(0, len(_s), 2): - _a, _b = _s.iloc[_k], _s.iloc[_k + 1] - if float(_a["p"]) > float(_b["p"]): # A = posizione con p minore - _a, _b = _b, _a - _righe.append({ - "pA": float(_a["p"]), "upA": float(_a["up"]), - "qA": float(_a["q"]), "uqA": float(_a["uq"]), - "pB": float(_b["p"]), "upB": float(_b["up"]), - "qB": float(_b["q"]), "uqB": float(_b["uq"]), - }) - _etichette.append(f"{_k}-{_k + 1}") - - coppie = pd.DataFrame(_righe, index=_etichette) - coppie - return (coppie,) - - -@app.cell -def _(Data, coppie, dati, sp): - _pA, _qA, _pB, _qB = sp.symbols("pA qA pB qB", positive=True) - _D = (_pA + _qA + _pB + _qB) / 2 - _d = (_pB - _pA + _qA - _qB) / 2 - - f3 = (_D**2 - _d**2) / (4 * _D) - - datiBessel = Data(dati.campione, analisi_stat=coppie.copy()) - datiBessel.calc_var(_D, "D") - datiBessel.calc_var(_d, "d") - datiBessel.calc_var(f3, "f3") - - res_3 = datiBessel.analisi_stat[ - ["D", "uD", "d", "ud", "f3", "uf3", "%uf3_pA", "%uf3_qA", "%uf3_pB", "%uf3_qB"] - ] - res_3.index = coppie.index - res_3 - return (res_3,) - - -@app.cell(hide_code=True) -def _(mo): - mo.md(r""" - ### Controlli di consistenza del metodo di Bessel - - - $D$ deve essere uguale nelle due righe: $p_A+q_A$ vs $p_B+q_B$ - - $d$ misurato con $p$ ($p_B-p_A$) e con $q$ ($q_A-q_B$) deve coincidere - - condizione di applicabilità $D>4f$ - - Compatibilità con `compat` (k < 2). - """) - return - - -@app.cell -def _(compat, coppie, np, pd, res_3): - _c = coppie - _DA, _uDA = _c.pA + _c.qA, np.sqrt(_c.upA**2 + _c.uqA**2) - _DB, _uDB = _c.pB + _c.qB, np.sqrt(_c.upB**2 + _c.uqB**2) - _dp, _udp = _c.pB - _c.pA, np.sqrt(_c.upA**2 + _c.upB**2) - _dq, _udq = _c.qA - _c.qB, np.sqrt(_c.uqA**2 + _c.uqB**2) - - _kD, _ = compat(_DA, _DB, _uDA, _uDB) - _kd, _ = compat(_dp, _dq, _udp, _udq) - - controlli_bessel = pd.DataFrame({ - "D_A": _DA, "D_B": _DB, "k(D)": _kD, "D compatibili": _kD < 2, - "d da p": _dp, "d da q": _dq, "k(d)": _kd, "d compatibili": _kd < 2, - "D > 4f": res_3["D"] > 4 * res_3["f3"], - }) - controlli_bessel - return - - -@app.cell(hide_code=True) -def _(mo): - mo.md(r""" - ## Risultati finali - - Per ogni metodo si combinano le stime (una per configurazione; una per - coppia nel metodo 3) con la **media pesata** ($w_i = 1/u_i^2$). Si - riportano anche $\chi^2$ e $P(\chi^2)$ di compatibilità delle stime con - un valore costante: se $\chi^2/\mathrm{dof}>1$ l'incertezza va ingrandita - di $\sqrt{\chi^2/\mathrm{dof}}$ (colonna `u_corr`). - """) - return - - -@app.cell -def _(np, pd, res_1, res_2a, res_2b, res_3, stats): - def media_pesata(f, uf): - f, uf = np.asarray(f, float), np.asarray(uf, float) - w = 1 / uf**2 - fm = np.sum(w * f) / np.sum(w) - ufm = 1 / np.sqrt(np.sum(w)) - chi2 = np.sum(w * (f - fm) ** 2) - dof = len(f) - 1 - return fm, ufm, chi2, dof, stats.chi2.sf(chi2, dof) - - _righe = {} - for _nome, _res, _c in [ - ("1 - punti coniugati", res_1, "f1"), - ("2a - p e ingrandimento", res_2a, "f2a"), - ("2b - q e ingrandimento", res_2b, "f2b"), - ("3 - Bessel", res_3, "f3"), - ]: - _fm, _u, _chi2, _dof, _P = media_pesata(_res[_c], _res[f"u{_c}"]) - _righe[_nome] = { - "f [cm]": _fm, "u": _u, - "chi2": _chi2, "dof": _dof, "P": _P, - "u_corr": _u * np.sqrt(max(1.0, _chi2 / _dof)), - } - - risultati = pd.DataFrame(_righe).T - risultati - return (risultati,) - - -@app.cell(hide_code=True) -def _(mo): - mo.md(r""" - ## Compatibilità tra i metodi e confronto grafico - """) - return - - -@app.cell -def _(compat, pd, risultati): - _nomi = list(risultati.index) - _righe = [] - for _i in range(len(_nomi)): - for _j in range(_i + 1, len(_nomi)): - _a, _b = risultati.loc[_nomi[_i]], risultati.loc[_nomi[_j]] - _k, _diff = compat(_a["f [cm]"], _b["f [cm]"], _a["u_corr"], _b["u_corr"]) - _righe.append({"metodi": f"{_nomi[_i][:2].strip()} vs {_nomi[_j][:2].strip()}", - "differenza": _diff, "k": _k, - "compatibili (k<2)": bool(_k < 2)}) - pd.DataFrame(_righe).set_index("metodi") - return - - -@app.cell -def _(plot_gauss, risultati): plot_gauss([ - (r["f [cm]"], r["u_corr"], _i, _n.strip()) - for _i, (_n, r) in enumerate(risultati.iterrows()) - ]) + ( + risultati_finali.loc["48B - m1", "f [cm]"], + risultati_finali.loc["48B - m1", "u_corr [cm]"], + 0, + "48B m1", + ), + ( + risultati_finali.loc["48B - m2", "f [cm]"], + risultati_finali.loc["48B - m2", "u_corr [cm]"], + 1, + "48B m2", + ), + ( + risultati_finali.loc["48B - m3", "f [cm]"], + risultati_finali.loc["48B - m3", "u_corr [cm]"], + 2, + "48B m3", + ) + ], cartella_immagini / "1.png", "Lente 48B - correzione angolare") + + plot_gauss([ + ( + risultati_finali.loc["252B - m1", "f [cm]"], + risultati_finali.loc["252B - m1", "u_corr [cm]"], + 3, + "252B m1", + ), + ( + risultati_finali.loc["252B - m2", "f [cm]"], + risultati_finali.loc["252B - m2", "u_corr [cm]"], + 4, + "252B m2", + ), + ( + risultati_finali.loc["252B - m3", "f [cm]"], + risultati_finali.loc["252B - m3", "u_corr [cm]"], + 5, + "252B m3", + ), + ], cartella_immagini / "2.png", "Lente 252B - correzione angolare") + + plot_gauss([ + ( + risultati_finali.loc["48B - m1", "f [cm]"], + risultati_finali.loc["48B - m1", "u_corr [cm]"], + 0, + "48B m1", + ), + (f48B_m2_senza, u48B_m2_senza, 1, "48B m2 senza correzione"), + ( + risultati_finali.loc["48B - m3", "f [cm]"], + risultati_finali.loc["48B - m3", "u_corr [cm]"], + 2, + "48B m3", + ), + ], cartella_immagini / "3.png", "Lente 48B - senza correzione angolare") + + plot_gauss([ + ( + risultati_finali.loc["252B - m1", "f [cm]"], + risultati_finali.loc["252B - m1", "u_corr [cm]"], + 3, + "252B m1", + ), + (f252B_m2_senza, u252B_m2_senza, 4, "252B m2 senza correzione"), + ( + risultati_finali.loc["252B - m3", "f [cm]"], + risultati_finali.loc["252B - m3", "u_corr [cm]"], + 5, + "252B m3", + ), + ], cartella_immagini / "4.png", "Lente 252B - senza correzione angolare") + return + + +@app.cell +def _(compat, pd, risultati_finali): + k_48B_m1_m2, diff_48B_m1_m2 = compat( + risultati_finali.loc["48B - m1", "f [cm]"], + risultati_finali.loc["48B - m2", "f [cm]"], + risultati_finali.loc["48B - m1", "u_corr [cm]"], + risultati_finali.loc["48B - m2", "u_corr [cm]"], + ) + k_48B_m1_m3, diff_48B_m1_m3 = compat( + risultati_finali.loc["48B - m1", "f [cm]"], + risultati_finali.loc["48B - m3", "f [cm]"], + risultati_finali.loc["48B - m1", "u_corr [cm]"], + risultati_finali.loc["48B - m3", "u_corr [cm]"], + ) + k_48B_m2_m3, diff_48B_m2_m3 = compat( + risultati_finali.loc["48B - m2", "f [cm]"], + risultati_finali.loc["48B - m3", "f [cm]"], + risultati_finali.loc["48B - m2", "u_corr [cm]"], + risultati_finali.loc["48B - m3", "u_corr [cm]"], + ) + k_252B_m1_m2, diff_252B_m1_m2 = compat( + risultati_finali.loc["252B - m1", "f [cm]"], + risultati_finali.loc["252B - m2", "f [cm]"], + risultati_finali.loc["252B - m1", "u_corr [cm]"] , + risultati_finali.loc["252B - m2", "u_corr [cm]"] , + ) + k_252B_m1_m3, diff_252B_m1_m3 = compat( + risultati_finali.loc["252B - m1", "f [cm]"], + risultati_finali.loc["252B - m3", "f [cm]"], + risultati_finali.loc["252B - m1", "u_corr [cm]"] , + risultati_finali.loc["252B - m3", "u_corr [cm]"] , + ) + k_252B_m2_m3, diff_252B_m2_m3 = compat( + risultati_finali.loc["252B - m2", "f [cm]"], + risultati_finali.loc["252B - m3", "f [cm]"], + risultati_finali.loc["252B - m2", "u_corr [cm]"] , + risultati_finali.loc["252B - m3", "u_corr [cm]"] , + ) + + confronti_metodi = pd.DataFrame({ + "48B m1 vs m2": [k_48B_m1_m2, diff_48B_m1_m2, k_48B_m1_m2 < 2], + "48B m1 vs m3": [k_48B_m1_m3, diff_48B_m1_m3, k_48B_m1_m3 < 2], + "48B m2 vs m3": [k_48B_m2_m3, diff_48B_m2_m3, k_48B_m2_m3 < 2], + "252B m1 vs m2": [k_252B_m1_m2, diff_252B_m1_m2, k_252B_m1_m2 < 2], + "252B m1 vs m3": [k_252B_m1_m3, diff_252B_m1_m3, k_252B_m1_m3 < 2], + "252B m2 vs m3": [k_252B_m2_m3, diff_252B_m2_m3, k_252B_m2_m3 < 2], + }, index=["k", "differenza [cm]", "compatibili (k < 2)"]).T + confronti_metodi + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + # Parte interessante + Ora ragioniamo sul perché avvengono certi fenomeni. + Risulta facile comprendere che il metodo 2 è quello più _sballato_ perché? + + Probabilmente abbiamo sofferto nel prendere h che subisce effetti di distorsione dovuti al fatto di avere lo schermo non fissato + una correzione di 15° è già sufficiente a sistemare il problema. Ovviamente oltre la correzione bisogna aggiungere degli errori su questa idea che sono $\sim$al valore stesso + + Sicuramente ci sono altri errori strumentali non considerati correttamente in quanto i $\chi^2$ son completamente fuori scala + """) return diff --git a/pyproject.toml b/pyproject.toml index 4989bfa..11c1fbb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,6 +5,7 @@ description = "Add your description here" readme = "README.md" requires-python = ">=3.14" dependencies = [ + "altair>=6.3.0", "marimo>=0.25.0", "matplotlib>=3.11.2", "numpy>=2.5.3", diff --git a/statistica.py b/statistica.py index 4ad75b5..17ac636 100644 --- a/statistica.py +++ b/statistica.py @@ -1,6 +1,6 @@ import marimo -__generated_with = "0.25.0" +__generated_with = "0.25.2" app = marimo.App() @@ -163,7 +163,7 @@ def _(np, pd, stats): return campione, rimossi - def stat_chauv(self, prefissi, err_strumentale): + def stat_chauv(self, prefissi, err_strumentale, err_della_media = True): campione_df = self.campione # riferimento al DataFrame originale (non sovrascritto) if isinstance(prefissi, str): @@ -185,15 +185,16 @@ def _(np, pd, stats): valori_puliti, rimossi = rimuovi_outlier(valori, err_strum) n = len(valori_puliti) - s = valori_puliti.std(ddof=1) / np.sqrt(n) + s = valori_puliti.std(ddof=1) + if err_della_media: s /= np.sqrt(n) u = np.sqrt(s**2 + err_strum**2) - self.analisi_stat.at[i, prefix] = valori_puliti.mean() - self.analisi_stat.at[i, f"u{prefix}_strum"] = err_strum - self.analisi_stat.at[i, f"u{prefix}_stat"] = s - self.analisi_stat.at[i, f"u{prefix}"] = u - self.analisi_stat.at[i, f"n{prefix}"] = n - self.analisi_stat.at[i, f"out{prefix}"] = rimossi + self.analisi_stat.loc[i, prefix] = valori_puliti.mean() + self.analisi_stat.loc[i, f"u{prefix}_strum"] = err_strum + self.analisi_stat.loc[i, f"u{prefix}_stat"] = s + self.analisi_stat.loc[i, f"u{prefix}"] = u + self.analisi_stat.loc[i, f"n{prefix}"] = n + self.analisi_stat.loc[i, f"out{prefix}"] = rimossi return self @@ -201,6 +202,95 @@ def _(np, pd, stats): return +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + # Media pesata + + input: + - array x: valori + - array ux: uvalori + output: + - xm: x medio + - uxm: sigma x medio + - chi2: $\chi^2$ + - dof: gradi di libertà + - Pchi2: $P(\chi^2)$ + + data: + - prefix: prefisso del nome della variabile + """) + return + + +@app.cell +def _(np, stats): + def media_pesata(x, ux): + x, ux = np.asarray(x, float), np.asarray(ux, float) + w = 1 / ux**2 + xm = np.sum(w * x) / np.sum(w) + uxm = 1 / np.sqrt(np.sum(w)) + chi2 = np.sum(w * (x - xm) ** 2) + dof = len(x) - 1 + return xm, uxm, chi2/dof, dof, stats.chi2.sf(chi2, dof) + + def media_pesata_data(self, prefix): + return media_pesata(self.analisi_stat[prefix], self.analisi_stat[f"u{prefix}"]) + + Data.media_pesata = media_pesata_data + return + + +@app.cell +def _(np, stats): + def media_pesata_dersimonian_laird(x, ux): + x, ux = np.asarray(x, float), np.asarray(ux, float) + n = len(x) + dof = n - 1 + + if n <= 1: + raise ValueError("Il dataset deve contenere almeno due misure.") + if np.any(ux <= 0): + raise ValueError("Le incertezze (ux) devono essere strettamente maggiori di zero.") + + # --- STEP 1: Calcolo standard (Effetti Fissi) per stimare il Chi2 originale --- + w_fixed = 1 / ux**2 + sum_w_fixed = np.sum(w_fixed) + xm_fixed = np.sum(w_fixed * x) / sum_w_fixed + chi2_fixed = np.sum(w_fixed * (x - xm_fixed) ** 2) + + # --- STEP 2: Stima della varianza extra (Tau^2) --- + # Formula di DerSimonian-Laird + sum_w2_fixed = np.sum(w_fixed**2) + numerator = chi2_fixed - dof + denominator = sum_w_fixed - (sum_w2_fixed / sum_w_fixed) + + # Tau^2 non può essere negativo. Se il Chi2 è minore del dof, Tau^2 = 0 + tau2 = max(0.0, numerator / denominator) + + # --- STEP 3: Ricalcolo con Effetti Casuali (Random Effects) --- + w_random = 1 / (ux**2 + tau2) + sum_w_random = np.sum(w_random) + + xm_random = np.sum(w_random * x) / sum_w_random + uxm_random = 1 / np.sqrt(sum_w_random) + + # Nota: Il p-value finale del Chi2 si riferisce alla consistenza iniziale dei dati + p_value = stats.chi2.sf(chi2_fixed, dof) + + return { + "media": xm_random, + "errore_media": uxm_random, + "tau2": tau2, + "tau": np.sqrt(tau2), + "chi2_originale": chi2_fixed, + "dof": dof, + "p_value": p_value + } + + return + + @app.cell(hide_code=True) def _(mo): mo.md(r""" @@ -589,7 +679,7 @@ def _(mo): @app.cell def _(np, plt, sns, stats): - def plot_gauss(gaussiane): + def plot_gauss(gaussiane, percorso=None, titolo="Confronto dei metodi"): # Creazione figura plt.figure(figsize=(12, 7)) @@ -600,10 +690,10 @@ def _(np, plt, sns, stats): for mu, sigma, _, _ in gaussiane: minimoLocale = mu - 4 * sigma massimoLocale = mu + 4 * sigma - + if minimoLocale < xMin: xMin = minimoLocale - + if massimoLocale > xMax: xMax = massimoLocale @@ -613,7 +703,7 @@ def _(np, plt, sns, stats): for mu, sigma, colore, etichetta in gaussiane: y = stats.norm.pdf(x, mu, sigma) plt.plot(x, y, color=sns.color_palette()[colore], linewidth=1, label=etichetta) - + puntiLinee = [mu - sigma, mu, mu + sigma] for px in puntiLinee: py = stats.norm.pdf(px, mu, sigma) @@ -621,13 +711,17 @@ def _(np, plt, sns, stats): ymin=0, ymax=py, colors=sns.color_palette()[colore], linestyles='dashed', linewidth=1) - + # Dettagli estetici finali plt.ylim(bottom=0) - plt.title('Confronto dati') - plt.xlabel('k') + plt.title(titolo) + plt.xlabel('Focale Lente [$cm$]') + plt.ylabel('Densità di probabilità [$cm^{-1}$]') plt.legend() + if percorso is not None: + plt.savefig(str(percorso), dpi=200, bbox_inches="tight") + # Mostriamo il grafico plt.show()