+1

analisis/plots/20231123_CPTconmicromocion3/lolo_analisis_superajuste.py

@@ -33,6 +33,87 @@ from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
 #     print(f'loaded: {var_name}')
 
 
+# Funciones auxiliares
+
+from scipy.stats.distributions import  t,chi2
+
+def estadistica(datos_x,datos_y,modelo,pcov,parametros,nombres=None,alpha=0.05):
+    
+    if nombres is None:
+        nombres = [ f'{j}' for j in range(len(parametros)) ]
+    
+    # Cantidad de parámetros
+    P = len(parametros)
+    
+    # Número de datos
+    N = len(datos_x)
+    
+    # Grados de libertas (Degrees Of Freedom)
+    dof = N-P-1
+
+    # Cauculamos coordenadas del modelo
+    # modelo_x    = datos_x if modelo_x_arr is None else modelo_x_arr
+    # modelo_y    = modelo( modelo_x, *parametros )
+    
+    # Predicción del modelo para los datos_x medidos
+    prediccion_modelo = modelo( datos_x, *parametros )
+    
+    # Calculos de cantidades estadísticas relevantes
+    COV       = pcov                                      # Matriz de Covarianza
+    SE        = np.sqrt(np.diag( COV  ))                        # Standar Error / Error estandar de los parámetros
+    residuos  = datos_y - prediccion_modelo               # diferencia enrte el modelo y los datos
+    
+    SSE       = sum(( residuos )**2 )                     # Resitual Sum of Squares
+    SST       = sum(( datos_y - np.mean(datos_y))**2)        # Total Sum of Squares
+    
+    # http://en.wikipedia.org/wiki/Coefficient_of_determination
+    # Expresa el porcentaje de la varianza que logra explicar el modelos propuesto
+    Rsq       =  1 - SSE/SST                               # Coeficiente de determinación
+    Rsq_adj   = 1-(1-Rsq) * (N-1)/(N-P-1)                  # Coeficiente de determinación Ajustado   
+    
+    # https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#In_least_squares_regression_analysis
+    # Expresa la correlación que hay entre los datos y la predicción del modelo
+    r_pearson = np.corrcoef( datos_y ,  prediccion_modelo )[0,1]
+    
+    # Reduced chi squared
+    # https://en.wikipedia.org/wiki/Reduced_chi-squared_statistic
+    chi2_     = sum( residuos**2 )/N
+    chi2_red  = sum( residuos**2 )/(N-P)
+    
+    # Chi squared test
+    chi2_test = sum( residuos**2 / abs(prediccion_modelo) )
+    # p-value del ajuste
+    p_val  = chi2(dof).cdf( chi2_test )
+    
+    
+    sT = t.ppf(1.0 - alpha/2.0, N - P ) # student T multiplier
+    CI = sT * SE                        # Confidence Interval
+    
+    print('R-squared    ',Rsq)
+    print('R-sq_adjusted',Rsq_adj)
+    print('chi2         ',chi2_)
+    print('chi2_reduced ',chi2_red)
+    print('chi2_test    ',chi2_test)
+    print('r-pearson    ',r_pearson)
+    print('p-value      ',p_val)
+    print('')
+    print('Error Estandard (SE):')
+    for i in range(P):
+        print(f'parametro[{nombres[i]:>5s}]: ' , parametros[i], ' ± ' , SE[i])
+    print('')
+    print('Intervalo de confianza al '+str((1-alpha)*100)+'%:')
+    for i in range(P):
+        print(f'parametro[{nombres[i]:>5s}]: ' , parametros[i], ' ± ' , CI[i])
+    
+    return dict(R2=Rsq,R2_adj=Rsq_adj,chi2=chi2_,chi2_red=chi2_red,
+                chi2_test=chi2_test,r=r_pearson,pvalue=p_val,
+                SE=SE,CI=CI)
+
+
+
+
+
+
 #%%
 
 """
@@ -160,6 +241,7 @@ def FitEIT_MM_single(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, T
 
 
 
+param_names = 'offset DetDoppler SG SP SCALE1 OFFSET BETA1 TEMP'.split()
 
 #%%
 """
@@ -233,6 +315,8 @@ if not 'popt_SA_vec' in globals().keys() or len(popt_SA_vec)==0:
 
     DetuningsUV_vec = []
     ErrorDetuningsUV_vec = []
+    
+    Estadistica_vec = []
 
     for selectedcurve in SelectedCurveVec:
         print(f"{round(time()-t0,1):6.1f}: Procesando la curva {selectedcurve}")
@@ -266,6 +350,12 @@ if not 'popt_SA_vec' in globals().keys() or len(popt_SA_vec)==0:
             FittedEITpi_3_SA_short, Detunings_3_SA_short = FitEIT_MM_single_plot(FreqsDR, *popt_3_SA)
             freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
             FittedEITpi_3_SA_long, Detunings_3_SA_long = FitEIT_MM_single_plot(freqslong, *popt_3_SA)
+            
+            
+            # estadistica(datos_x,datos_y,modelo,pcov,parametros,nombres=None,alpha=0.05)
+            est_tmp = estadistica(FreqsDR,CountsDR,FitEIT_MM_single,pcov_3_SA,popt_3_SA,
+                                  nombres=param_names,alpha=0.05)
+            Estadistica_vec.append(est_tmp)
 
         DetuningsUV_vec.append(popt_3_SA[1])
         ErrorDetuningsUV_vec.append(np.sqrt(pcov_3_SA[1,1]))
@@ -335,7 +425,10 @@ param_names = 'offset DetDoppler SG SP SCALE1 OFFSET BETA1 TEMP'.split()
 err_vecs = np.array([ np.sqrt(np.diag(el)) for el in pcov_SA_vec ])
 num_med  = np.arange(len(pcov_SA_vec)) +1
 
-fig, axx = plt.subplots( len(popt_SA_vec[0]),1, figsize=(13,8) ,  constrained_layout=True, sharex=True , sharey=False )
+r2_values = np.array([ el['R2_adj'] for el in Estadistica_vec ])
+
+
+fig, axx = plt.subplots( len(popt_SA_vec[0])+1,1, figsize=(13,8) ,  constrained_layout=True, sharex=True , sharey=False )
 fig.set_constrained_layout_pads(w_pad=2/72, h_pad=2/72, hspace=0, wspace=0)
 
 for ax,param_vec,err_vec,par_name in zip(axx,popt_SA_vec.T,err_vecs.T,param_names) :
@@ -346,6 +439,15 @@ for ax,param_vec,err_vec,par_name in zip(axx,popt_SA_vec.T,err_vecs.T,param_name
     ax.grid(True, ls=":", color='lightgray')
     ax.set_ylabel(par_name)
 
+
+ax=axx[-1]
+ax.plot( num_med , r2_values, '.-')
+ax.set_ylabel(r'$R^2$')
+ax.grid(True, ls=":", color='lightgray')
+
+
+
+fig.align_ylabels()
 ax.set_xticks(num_med)
 ax.set_xlabel('Num. de medición')