prueba_t_para_una_muestra.py

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from typing import Tuple, Callable

class PowerFunctionOneSampleT:
    """
    Calcula y visualiza la función de poder para prueba t de una muestra.
    
    Conexión con Sesión 1:
    - Espacio muestral: Ω = ℝⁿ con σ-álgebra de Borel ℬ(ℝⁿ)
    - Medida: ℙ_μ = 𝒩(μ, σ²)⊗ⁿ (medida producto)
    - Regla φ: función medible φ: ℝⁿ → {0,1}
    """
    
    def __init__(self, n: int, sigma: float, mu_0: float, alpha: float = 0.05):
        self.n = n
        self.sigma = sigma
        self.mu_0 = mu_0
        self.alpha = alpha
        self.df = n - 1  # grados de libertad
        self.t_critical = stats.t.ppf(1 - alpha/2, self.df)
        
    def noncentrality_parameter(self, mu: float) -> float:
        """Calcula el parámetro de no centralidad δ para un μ dado"""
        return (mu - self.mu_0) / (self.sigma / np.sqrt(self.n))
    
    def power_at_mu(self, mu: float) -> float:
        """
        Calcula β_φ(μ) = P_μ(Rechazar H₀)
        
        Fundamento: Bajo H₁, T ~ t_{n-1}(δ) con δ = (μ - μ₀)/(σ/√n)
        """
        delta = self.noncentrality_parameter(mu)
        
        # Potencia = P(T < -t_crit) + P(T > t_crit) bajo t no central
        power_lower = stats.nct.cdf(-self.t_critical, self.df, delta)
        power_upper = stats.nct.sf(self.t_critical, self.df, delta)
        
        return power_lower + power_upper
    
    def power_function(self, mu_values: np.ndarray) -> np.ndarray:
        """Evalúa β_φ(μ) para un array de valores de μ"""
        return np.array([self.power_at_mu(mu) for mu in mu_values])
    
    def type2_error_function(self, mu_values: np.ndarray) -> np.ndarray:
        """Calcula β(μ) = 1 - β_φ(μ) (Error Tipo II)"""
        return 1 - self.power_function(mu_values)
    
    def plot_power_curve(self, mu_range: Tuple[float, float], 
                         n_points: int = 500,
                         ax: plt.Axes = None) -> plt.Axes:
        """
        Visualiza la función de poder β_φ(μ)
        
        Elementos gráficos:
        - Curva de poder principal
        - Línea de nivel α en Θ₀
        - Región de potencia mínima (1-β) recomendada
        """
        if ax is None:
            fig, ax = plt.subplots(figsize=(12, 7))
        
        mu_values = np.linspace(mu_range[0], mu_range[1], n_points)
        power_values = self.power_function(mu_values)
        
        # Curva de poder
        ax.plot(mu_values, power_values, 'b-', linewidth=2.5, 
                label=r'Función de Poder $\beta_\phi(\mu)$')
        
        # Línea de nivel α (en H₀: μ = μ₀)
        ax.axvline(x=self.mu_0, color='red', linestyle='--', linewidth=2, 
                   label=r'$H_0: \mu = \mu_0$')
        ax.axhline(y=self.alpha, color='green', linestyle=':', linewidth=2,
                   label=r'Nivel $\alpha = {:.2f}$'.format(self.alpha))
        
        # Región de potencia recomendada (≥ 0.80)
        ax.axhline(y=0.80, color='orange', linestyle='-.', linewidth=2,
                   label='Potencia mínima recomendada (80%)')
        ax.fill_between(mu_values, 0.80, 1.0, alpha=0.2, color='orange',
                       label='Región de potencia adecuada')
        
        # Punto de potencia en μ₀ (debe ser ≈ α)
        power_at_null = self.power_at_mu(self.mu_0)
        ax.plot(self.mu_0, power_at_null, 'ro', markersize=10,
               label=r'$\beta_\phi(\mu_0) = {:.4f}$'.format(power_at_null))
        
        # Etiquetas y formato
        ax.set_xlabel(r'Parámetro $\mu$', fontsize=14)
        ax.set_ylabel(r'Probabilidad $\beta_\phi(\mu)$', fontsize=14)
        ax.set_title(r'Función de Poder para Prueba $t$ (n={}, $\sigma={:.2f}$)'.format(
            self.n, self.sigma), fontsize=14, fontweight='bold')
        ax.legend(loc='upper left', fontsize=11)
        ax.grid(alpha=0.3, linestyle='--')
        ax.set_ylim(0, 1.05)
        
        # Anotación del efecto mínimo detectable
        # Encontrar μ donde potencia = 0.80
        from scipy.optimize import brentq
        try:
            mu_80_right = brentq(lambda m: self.power_at_mu(m) - 0.80, 
                                self.mu_0, mu_range[1])
            mu_80_left = brentq(lambda m: self.power_at_mu(m) - 0.80, 
                               mu_range[0], self.mu_0)
            
            ax.annotate(f'MED: μ = {mu_80_right:.2f}', 
                       xy=(mu_80_right, 0.80), xytext=(mu_80_right+0.3, 0.75),
                       arrowprops=dict(arrowstyle='->', color='black'),
                       fontsize=10, fontweight='bold')
        except:
            pass  # Si no encuentra raíz, omitir anotación
        
        return ax
    
    def summary_table(self, mu_alternatives: np.ndarray) -> None:
        """Imprime tabla resumen de poder para valores específicos de μ"""
        print(f"{'μ':>10} | {'δ (ncp)':>10} | {'Potencia β_φ(μ)':>15} | {'Error Tipo II β(μ)':>15}")
        print("-" * 60)
        for mu in mu_alternatives:
            delta = self.noncentrality_parameter(mu)
            power = self.power_at_mu(mu)
            type2 = 1 - power
            print(f"{mu:>10.3f} | {delta:>10.3f} | {power:>15.4f} | {type2:>15.4f}")