Appendice D — Probabilità

# Standard library imports
import os

# Third-party imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import arviz as az
import scipy.stats as stats
from scipy.special import expit  # Funzione logistica
import math
from cmdstanpy import cmdstan_path, CmdStanModel

# Configuration
seed = sum(map(ord, "stan_poisson_regression"))
rng = np.random.default_rng(seed=seed)
az.style.use("arviz-darkgrid")
%config InlineBackend.figure_format = "retina"

# Define directories
home_directory = os.path.expanduser("~")
project_directory = f"{home_directory}/_repositories/psicometria"

# Print project directory to verify
print(f"Project directory: {project_directory}")

Project directory: /Users/corradocaudek/_repositories/psicometria

?sec-prob-on-general-spaces

?exr-prob-on-general-spaces-1

Per calcolare questa probabilità in maniera analitica, utilizziamo la seguente uguaglianza:

\[ P(\text{almeno 2 psicologi clinici}) = 1 - P(\text{nessun psicologo clinico}) - P(\text{1 psicologo clinico}). \]

Il numero totale di modi per selezionare 5 persone dal gruppo di 20 è dato da:

\[ \binom{20}{5} = \frac{20!}{5!(15!)} = 15,504. \]

Il numero di modi per avere nessun psicologo clinico nella commissione (ovvero, selezionare solo psicologi del lavoro) è:

\[ \binom{10}{0} \times \binom{10}{5} = 1 \times 252 = 252. \]

Quindi, la probabilità di avere nessun psicologo clinico è:

\[ P(\text{nessun psicologo clinico}) = \frac{252}{15,504} \approx 0.016. \]

Il numero di modi per avere esattamente 1 psicologo clinico nella commissione è:

\[ \binom{10}{1} \times \binom{10}{4} = 10 \times 210 = 2,100. \]

Quindi, la probabilità di avere esattamente 1 psicologo clinico è:

\[ P(\text{1 psicologo clinico}) = \frac{2,100}{15,504} \approx 0.135. \]

La probabilità di avere almeno 2 psicologi clinici nella commissione è quindi:

\[ \begin{align} P(\text{almeno 2 psicologi clinici}) &= 1 - P(\text{nessun psicologo clinico}) - P(\text{1 psicologo clinico}) \notag\\ &= 1 - 0.016 - 0.135 \notag\\ &= 0.848.\notag \end{align} \]

Quindi, la probabilità che almeno 2 psicologi clinici siano nella commissione è circa 0.848.

# Funzione per calcolare le combinazioni
def nCk(n, k):
    return math.factorial(n) // (math.factorial(k) * math.factorial(n - k))


# Calcolo delle probabilità per il problema della commissione
total_ways = nCk(20, 5)
no_clinical = nCk(10, 0) * nCk(10, 5)
one_clinical = nCk(10, 1) * nCk(10, 4)

p_no_clinical = no_clinical / total_ways
p_one_clinical = one_clinical / total_ways

p_at_least_two_clinical = 1 - p_no_clinical - p_one_clinical

print(f"Probabilità di almeno 2 psicologi clinici: {p_at_least_two_clinical:.3f}")

Probabilità di almeno 2 psicologi clinici: 0.848

In maniera più intuitiva, possiamo risolvere il problema con una simulazione Monte Carlo.

import random

# Numero di simulazioni
simulations = 1000000

# Numero di successi (almeno 2 psicologi clinici nella commissione)
success_count = 0

# Creiamo una lista che rappresenta il gruppo di 20 persone
# 1 rappresenta un psicologo clinico, 0 rappresenta un psicologo del lavoro
group = [1] * 10 + [0] * 10

# Simulazione Monte Carlo
for _ in range(simulations):
    # Estrai casualmente 5 persone dal gruppo
    committee = random.sample(group, 5)

    # Conta quanti psicologi clinici ci sono nella commissione
    num_clinical_psychologists = sum(committee)

    # Verifica se ci sono almeno 2 psicologi clinici
    if num_clinical_psychologists >= 2:
        success_count += 1

# Calcola la probabilità
probability = success_count / simulations

# Mostra il risultato
print(
    f"La probabilità che almeno 2 psicologi clinici siano nella commissione è: {probability:.4f}"
)

La probabilità che almeno 2 psicologi clinici siano nella commissione è: 0.8482

?sec-simulations

?exr-prob-simulation-1

Per calcolare le deviazioni standard delle distribuzioni gaussiane date le percentuali di studenti che ottengono meno di 18, possiamo utilizzare le proprietà della distribuzione normale e i quantili della distribuzione normale standard (distribuzione normale con media 0 e deviazione standard 1).

Le distribuzioni normali hanno la proprietà che possiamo trasformare qualsiasi valore \(X\) della distribuzione \(N(\mu, \sigma)\) nella distribuzione normale standard \(N(0, 1)\) tramite la formula:

\[ Z = \frac{X - \mu}{\sigma}, \]

dove \(Z\) è il quantile standardizzato.

Per trovare il valore di \(\sigma\) dato un certo percentile, utilizziamo l’inverso della funzione di distribuzione cumulativa (CDF) della distribuzione normale standard. Per un dato percentile \(p\), \(z_p\) è tale che:

\[ p = P(Z \leq z_p) \]

Quindi possiamo trovare \(\sigma\) risolvendo per \(\sigma\) nella formula:

\[ z_p = \frac{X - \mu}{\sigma}, \]

\[ \sigma = \frac{X - \mu}{z_p}, \]

dove:

\(X\) è il punteggio di soglia (18 in questo caso).
\(\mu\) è la media della distribuzione.
\(z_p\) è il quantile della distribuzione normale standard per il percentile \(p\).

I quantili della distribuzione normale standard per i percentili desiderati sono:

Per il 15%, il quantile è \(z_{0.15} \approx -1.036\).
Per il 10%, il quantile è \(z_{0.10} \approx -1.281\).
Per il 5%, il quantile è \(z_{0.05} \approx -1.645\).

Prima Prova

Media: \(\mu = 24\)
Percentuale che ottiene meno di 18: 15%
Quantile: \(z_{0.15} = -1.036\)
Soglia: \(X = 18\)

\[ \sigma_1 = \frac{24 - 18}{1.036} \approx 5.79 \]

Seconda Prova

Media: \(\mu = 25\)
Percentuale che ottiene meno di 18: 10%
Quantile: \(z_{0.10} = -1.281\)
Soglia: \(X = 18\)

\[ \sigma_2 = \frac{25 - 18}{1.281} \approx 5.46 \]

Terza Prova

Media: \(\mu = 26\)
Percentuale che ottiene meno di 18: 5%
Quantile: \(z_{0.05} = -1.645\)
Soglia: \(X = 18\)

\[ \sigma_3 = \frac{26 - 18}{1.645} \approx 4.86 \]

# Funzione per calcolare la deviazione standard data la media, la soglia e il quantile
def calculate_std(mean, threshold, quantile):
    return abs((mean - threshold) / quantile)


# Parametri delle distribuzioni gaussiane per le tre prove
mean_test1 = 24
std_test1 = calculate_std(mean_test1, 18, -1.036)
mean_test2 = 25
std_test2 = calculate_std(mean_test2, 18, -1.281)
mean_test3 = 26
std_test3 = calculate_std(mean_test3, 18, -1.645)

# Numero di studenti
n_students = 220

# Percentuale di studenti che non fa le prove
drop_test1 = 0.10
drop_test2 = 0.05

# Seed per il generatore di numeri casuali basato sulla stringa "simulation"
seed = sum(map(ord, "simulation"))
rng = np.random.default_rng(seed=seed)

# Generazione dei voti per le tre prove
# Genera i voti solo per gli studenti che partecipano alla prova
test1_scores = np.where(
    rng.random(n_students) > drop_test1,
    rng.normal(mean_test1, std_test1, n_students),
    np.nan,
)
test2_scores = np.where(
    rng.random(n_students) > drop_test2,
    rng.normal(mean_test2, std_test2, n_students),
    np.nan,
)
test3_scores = rng.normal(mean_test3, std_test3, n_students)

# Calcola il voto finale solo per gli studenti che hanno partecipato a tutte e tre le prove
final_scores = np.nanmean(
    np.column_stack((test1_scores, test2_scores, test3_scores)), axis=1
)

# Filtra gli studenti che non hanno partecipato a tutte e tre le prove
valid_final_scores = final_scores[~np.isnan(final_scores)]

# Visualizzazione della distribuzione finale dei voti
plt.hist(valid_final_scores, bins=30, edgecolor="black")
plt.title("Distribuzione dei voti finali")
plt.xlabel("Voto finale")
plt.ylabel("Frequenza")
plt.show()

# Statistiche descrittive dei voti finali
mean_final_score = np.mean(valid_final_scores)
median_final_score = np.median(valid_final_scores)
std_final_score = np.std(valid_final_scores)

print(f"Media dei voti finali: {mean_final_score:.2f}")
print(f"Mediana dei voti finali: {median_final_score:.2f}")
print(f"Deviazione standard dei voti finali: {std_final_score:.2f}")