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

# Definir la función
def trig_function(x, A, B, C):
    return A * np.sin(B * x + C)
# Definir la derivada de la función
def dfdx(x, A, B, C):
    return A * B * np.cos(B * x + C)

# Defino los parametros
Amplitud = 1
Frecuencia = 1
Fase = 0
peso = 1.5 # asumiendo que los errores en X fueron un 50% mas grandes

# Generar puntos en el rango de interes
x_values = np.linspace(0, 10*np.pi, 50)

# Calcular los valores que toma la funcion
y_true = trig_function(x_values, Amplitud, Frecuencia, Fase)

# Calcular los valores de la derivada de la función 
dx = dfdx(x_values, Amplitud, Frecuencia, Fase)

# Introducir errores en los puntos
error_porcentual_X = 0.20
error_porcentual_Y = 0.05 

error_x = np.random.normal(0, error_porcentual_X, len(x_values))
error_y = np.random.normal(0, error_porcentual_Y, len(x_values))

# Aplicar errores a los puntos
x_values_with_errors = x_values + error_x
y_values_with_errors = y_true + error_y

# Ajustar la función a los datos con errores en ambos ejes usando con curve_fit
popt, pcov = curve_fit(trig_function, x_values_with_errors, y_values_with_errors, absolute_sigma=True)

# Calcular los residuos
residuos = y_values_with_errors - trig_function(x_values_with_errors, *popt)

# Obtener los errores estándar de los parámetros ajustados
perr = np.sqrt(np.diag(pcov))

# Graficar los puntos con barras de error en x e y
plt.figure(figsize=(10, 6))
plt.errorbar(x_values_with_errors, y_values_with_errors, xerr=np.abs(error_x), yerr=np.abs(error_y), fmt='o', color='red', alpha=0.5, label='Puntos con Errores')
plt.plot(x_values, trig_function(x_values, *popt), label='Ajuste de la Función', color='green', linestyle='--')
plt.title('Puntos con Errores y Ajuste de la Función')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.grid(True)
plt.show()

# Calcular el chi-cuadrado
sigma_cuadrado = (peso *dx*error_x)**2 + (error_y**2)
chi2_val = np.sum(residuos**2 / sigma_cuadrado)

# Calcular los grados de libertad
df = len(x_values_with_errors) - len(popt)

# Calcular el p-valor
p_value = 1 - chi2.cdf(chi2_val, df)

# Imprimir resultados
print(f"Chi-cuadrado: {chi2_val:.4f}")
print(f"Grados de Libertad: {df:.4f}")
print(f"Valor p asociado: {p_value:.4f}")

print("Parámetros del ajuste:")
print(f"  A = {popt[0]:.4f} +/- {perr[0]:.4f}")
print(f"  B = {popt[1]:.4f} +/- {perr[1]:.4f}")
print(f"  C = {popt[2]:.4f} +/- {perr[2]:.4f}")