import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import chi2
import matplotlib.pyplot as plt

# Definimos la función sinusoidal
def sinusoidal(x, A, omega, phi, offset):
    return A * np.sin(omega * x + phi) + offset

# Generamos datos sin ruido, esto sera luego reemplazado por las mediciones
par_iniciales = [2.0, 1.5, 0.5, 1.0]  
np.random.seed(0)  # Fijamos la semilla para reproducibilidad
x_data = np.linspace(0, 10*np.pi, 500)  # simulamos 500 puntos
y_true = sinusoidal(x_data, par_iniciales[0], par_iniciales[1], par_iniciales[2], par_iniciales[3])

# Agregamos error a  las mediciones simuladas
y_data = y_true + np.random.normal(scale=0.1, size=x_data.size)

# Ajuste por mínimos cuadrados con parámetros iniciales estimados
# Valores iniciales sugeridos para el ajuste, los mismos con que simulamos los datos. 
# Si tenemos una funcion medida, estimar mirando los datos valores razonables aproximados para estos parametros.

popt, pcov = curve_fit(sinusoidal, x_data, y_data, p0=par_iniciales)

# Parámetros del ajuste y sus errores
A_fit, omega_fit, phi_fit, offset_fit = popt
A_error, omega_error, phi_error, offset_error = np.sqrt(np.diag(pcov))

# Cálculo del chi-cuadrado
y_fit = sinusoidal(x_data, *popt)
residuos = y_data - y_fit
sigma2=y_fit  # esto equivale a asumir que los errores son Poissonianos
chi_squared = np.sum(residuos**2 / sigma2)
degrees_of_freedom = len(x_data) - len(popt)
p_value = 1 - chi2.cdf(chi_squared, degrees_of_freedom)

# Reporte de resultados
print("Parámetros del ajuste:")
print("Amplitud (A): {:.3f} ± {:.3f}".format(A_fit, A_error))
print("Frecuencia (omega): {:.3f} ± {:.3f}".format(omega_fit, omega_error))
print("Fase (phi): {:.3f} ± {:.3f}".format(phi_fit, phi_error))
print("Offset: {:.3f} ± {:.3f}".format(offset_fit, offset_error))
print("\nChi-cuadrado:", chi_squared)
print("Grados de libertad:", degrees_of_freedom)
print("Valor-p:", p_value)

# Graficar datos y ajuste
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label='Datos con ruido')
plt.plot(x_data, y_true, 'g--', label='Función sinusoidal original')
plt.plot(x_data, y_fit, 'r-', label='Ajuste sinusoidal')
plt.xlabel('Tiempo [s]', fontsize=18)
plt.ylabel('Amplitud', fontsize=18)
plt.legend(fontsize=18)
plt.grid(True)
plt.show()