fix

4dc6556b · Tatiana Acero Cuellar · 2fe8629f · 4dc6556b · 4dc6556b · 4dc6556b
Commit 4dc6556b authored 4 years ago by Tatiana Acero Cuellar
--- a/Functions/functions.py
+++ b/Functions/functions.py
@@ -304,8 +304,8 @@ def grafica_data_ybestfit(path, gal, df_best_parmeters, ss="ra"):

    ax2 = plt.subplot(grilla[2:,:])
    # -- graficar residuo
-    ax2.plot(gal["r(kpc)"],(ynfw(gal["r(kpc)"]*1000) - gal["v(km/s)"]) / np.sqrt(gal["err_v(km/s)"]) ,'--r',color="darkorange")
-    ax2.plot(gal["r(kpc)"],(yiso(gal["r(kpc)"]*1000) - gal["v(km/s)"]) / np.sqrt(gal["err_v(km/s)"]) ,'--r',color="blue")
+    ax2.plot(gal["r(kpc)"],(ynfw(gal["r(kpc)"]*1000) - gal["v(km/s)"]) / (np.sqrt(gal["err_v(km/s)"])**2) ,'--r',color="darkorange")
+    ax2.plot(gal["r(kpc)"],(yiso(gal["r(kpc)"]*1000) - gal["v(km/s)"]) / (np.sqrt(gal["err_v(km/s)"])**2) ,'--r',color="blue")

        
    # -- poner nombre a ejes y título a la gráfica

--- a/notebook/confidence_interval_ellipse1-short - 1.ipynb
+++ b/notebook/confidence_interval_ellipse1-short - 1.ipynb
@@ -278,6 +278,118 @@
    "ch2_inter_ISO_ = chi2_intervals(best_parmeters, df1, chis_ISO_ra, i = 1, ss = \"ra\")"
   ]
  },
+  {
+   "cell_type": "code",
+   "execution_count": 56,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m01 = np.exp(-np.array(list(zip(*ch2_inter_NFW_[1]))[1])/2).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 57,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m12 = np.exp(-np.array(list(zip(*ch2_inter_NFW_[2]))[1])/2).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 58,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m23 = np.exp(-np.array(list(zip(*ch2_inter_NFW_[3]))[1])/2).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 59,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m33 = np.exp(-np.array(list(zip(*ch2_inter_NFW_[4]))[1])/2).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 61,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "total = m01+m12+m23+m33"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 69,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.09570840914791515"
+      ]
+     },
+     "execution_count": 69,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "m01 / total"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 63,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.16337829265551743"
+      ]
+     },
+     "execution_count": 63,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "m12 / total"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 67,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "6.147853929639663e-97"
+      ]
+     },
+     "execution_count": 67,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.exp(-chi2_NFW_bestfit/2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
  {
   "cell_type": "code",
   "execution_count": 10,
@@ -382,16 +494,6 @@
    "plot_ellipse_confidence_interval(path,para_range_ISO, best_parmeters, ss=\"ra\", i=1 )"
   ]
  },
-  {
-   "cell_type": "code",
-   "execution_count": 16,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# -- gráfica, en functions.py, toca descomentar la línea que guarda la gráfica\n",
-    "#grafica_data_ybestfit(path, df1, best_parmeters, ss=\"ra\", mtype = \"new\", new_para_min_ch2_b = new_para_min_ch2_b)"
-   ]
-  },
  {
   "cell_type": "code",
   "execution_count": 17,
@@ -407,10 +509,10 @@
    "# ,len(df1) # saca la longitud de los datos\n",
    "# ,best_parmeters[\"M200_NFW(Msun)\"][0] # aca toca cambiarle el 0 (ra), por 1(r) o 2(a)\n",
    "# ,best_parmeters[\"c_NFW\"][0] # aca toca cambiarle el 0 (ra), por 1(r) o 2(a)\n",
-    "# ,chi2_NFW_bestfit # chi2 para NFW con los parámetros que da leastsq\n",
+    "# ,chi2_NFW_bestfit_red # chi2 para NFW con los parámetros que da leastsq\n",
    "# ,abs(best_parmeters[\"rho0_ISO(Msun/pc3)\"][0]) # aaca toca cambiarle el 0 (ra), por 1(r) o 2(a)\n",
    "# ,abs(best_parmeters[\"Rc_ISO(pc)\"][0]) # aca toca cambiarle el 0 (ra), por 1(r) o 2(a)\n",
-    "# ,chi2_ISO_bestfit)) # chi2 para ISO con los parámetros que da leastsq\n",
+    "# ,chi2_ISO_bestfit_red)) # chi2 para ISO con los parámetros que da leastsq\n",
    "# file1.close()\n"
   ]
  },

 %% Cell type:code id: tags:

 ``` python
 import sys
 # insert at 1, 0 is the script path (or '' in REPL)
 sys.path.insert(1, '../Functions')

 from imports import *
 from constants import *
 from functions import *
 ```

 %% Cell type:code id: tags:

 ``` python
 # path = ["../galaxies_data/rcugc4165.txt"]
 path = "../galaxies_data/rcugc4165.txt"
 ```

 %% Cell type:code id: tags:

 ``` python
 AA = []
 #for p in path:
 print(path)
 col = read_col_names(path)
 df1 = pd.read_csv(path,skiprows=1,header=None,delimiter=r"\s+",names=col)
 # -- quitar errores nulos para evitar problemas con leastsq
 df1 = df1[df1["err_v(km/s)"]!=0]
 print(len(df1))

 for labe in ['ra','r','a']:
    if labe == 'ra':
        gal = df1.copy()
    if labe == 'r':
        dfr = df1.copy()[df1.copy()['Side'] == 'r']
        gal = dfr
    if labe == 'a':
        dfa = df1.copy()[df1.copy()['Side'] == 'a']
        gal = dfa

    AA.append((path.split("/")[-1][:-4],determine(gal,labe)))
 best_parmeters = return_values(AA)
 best_parmeters
 ```

 %% Output

    ../galaxies_data/rcugc4165.txt
    128

    ../Functions/functions.py:67: RuntimeWarning: invalid value encountered in sqrt
      v = ( 4*np.pi*G*rr*x**(-1) *(x - np.arctan(x) ) )**(1/2)
    ../Functions/functions.py:202: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
      A = np.array(AA)

              ID Side  M200_NFW(Msun)  err_M200_NFW(Msun)      c_NFW  err_c_NFW  \
    0  rcugc4165   ra    6.074409e+10        1.228760e+10  19.137209   1.851499
    1  rcugc4165    r    3.300000e+10        6.069498e+09  28.515061   3.587851
    2  rcugc4165    a    3.562725e+11        2.543976e+11   9.808096   2.170091
    
       rho0_ISO(Msun/pc3)  err_rho0_ISO(Msun/pc3)  Rc_ISO(pc)  err_Rc_ISO(pc)
    0            0.394941                0.073572  564.171234       65.353228
    1            0.979249                0.371557 -341.529478       72.642947
    2            0.174522                0.027822 -967.563968      118.061555

 %% Cell type:code id: tags:

 ``` python
 # -- gráfica, en functions.py, toca descomentar la línea que guarda la gráfica
 grafica_data_ybestfit(path, df1, best_parmeters, ss="ra")
 ```

 %% Output



 %% Cell type:code id: tags:

 ``` python
 # aca toca cambiarle el 0 (ra), por 1(r) o 2(a)
 y_pred_NFW_bestfit = vnfw(best_parmeters["M200_NFW(Msun)"][0], best_parmeters["c_NFW"][0], df1["r(kpc)"]*1000/r200(best_parmeters["M200_NFW(Msun)"][0]), G = G, H0 = H0)
 y_pred_ISO_bestfit = viso((best_parmeters["rho0_ISO(Msun/pc3)"][0]*best_parmeters["Rc_ISO(pc)"][0]**2), df1["r(kpc)"]*1000/abs(best_parmeters["Rc_ISO(pc)"][0]), G = G)


 chi2_NFW_bestfit = (((y_pred_NFW_bestfit-df1["v(km/s)"])/(df1["err_v(km/s)"]))**2).sum()
 chi2_ISO_bestfit = (((y_pred_ISO_bestfit-df1["v(km/s)"])/(df1["err_v(km/s)"]))**2).sum()
 print("Mejor $\chi^2$ para \nNFW: {} \nISO: {}".format(chi2_NFW_bestfit,chi2_ISO_bestfit))

 chi2_NFW_bestfit_red = chi2_NFW_bestfit / (len(df1)-2)
 chi2_ISO_bestfit_red = chi2_ISO_bestfit / (len(df1)-2)
 print("Mejor $\chi^2$ reducido para \nNFW: {} \nISO: {}".format(chi2_NFW_bestfit_red,chi2_ISO_bestfit_red))

 ```

 %% Output

    Mejor $\chi^2$ para
    NFW: 443.0693019080659
    ISO: 464.11643927074226
    Mejor $\chi^2$ reducido para
    NFW: 3.5164230310163958
    ISO: 3.68346380373605

 %% Cell type:code id: tags:

 ``` python
 print("La probabilidad cumulativa de encontrar un $\chi^2$ menor para NFW es: {}".format(scipy.stats.chi2.cdf(chi2_NFW_bestfit,len(df1)-2)))
 print("La probabilidad cumulativa de encontrar un $\chi^2$ menor para ISO es: {}".format(scipy.stats.chi2.cdf(chi2_ISO_bestfit,len(df1)-2)))

 ```

 %% Output

    La probabilidad cumulativa de encontrar un $\chi^2$ menor para NFW es: 1.0
    La probabilidad cumulativa de encontrar un $\chi^2$ menor para ISO es: 1.0

 %% Cell type:code id: tags:

 ``` python
 # -- i = 0: NFW
 combs_NFW_ra = parameters_range(best_parmeters, i = 0, ss = "ra")
 # -- i = 1: ISO
 combs_ISO_ra = parameters_range(best_parmeters, i = 1, ss = "ra")
 ```

 %% Cell type:code id: tags:

 ``` python
 # -- i = 0: NFW
 chis_NFW_ra = chi2_parameters_range(df1, combs_NFW_ra, i = 0)
 # -- i = 1: ISO
 chis_ISO_ra = chi2_parameters_range(df1, combs_ISO_ra, i = 1)
 ```

 %% Cell type:code id: tags:

 ``` python
 ch2_inter_NFW_  = chi2_intervals(best_parmeters, df1, chis_NFW_ra, i = 0, ss = "ra")
 ch2_inter_ISO_ = chi2_intervals(best_parmeters, df1, chis_ISO_ra, i = 1, ss = "ra")
 ```

 %% Cell type:code id: tags:

 ``` python
+m01 = np.exp(-np.array(list(zip(*ch2_inter_NFW_[1]))[1])/2).sum()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+m12 = np.exp(-np.array(list(zip(*ch2_inter_NFW_[2]))[1])/2).sum()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+m23 = np.exp(-np.array(list(zip(*ch2_inter_NFW_[3]))[1])/2).sum()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+m33 = np.exp(-np.array(list(zip(*ch2_inter_NFW_[4]))[1])/2).sum()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+total = m01+m12+m23+m33
+```
+
+%% Cell type:code id: tags:
+
+``` python
+m01 / total
+```
+
+%% Output
+
+    0.09570840914791515
+
+%% Cell type:code id: tags:
+
+``` python
+m12 / total
+```
+
+%% Output
+
+    0.16337829265551743
+
+%% Cell type:code id: tags:
+
+``` python
+np.exp(-chi2_NFW_bestfit/2)
+```
+
+%% Output
+
+    6.147853929639663e-97
+
+%% Cell type:code id: tags:
+
+``` python
+```
+
+%% Cell type:code id: tags:
+
+``` python
 para_range_NFW = []
 para_range_ISO = []
 for kkk in range(1,5):
    para_range_NFW.append(parameters_each_interval(combs_NFW_ra, ch2_inter_NFW_, kk = kkk))
    para_range_ISO.append(parameters_each_interval(combs_ISO_ra, ch2_inter_ISO_, kk = kkk))
 ```

 %% Cell type:code id: tags:

 ``` python
 # -- índice dentro del arreglo de chi^2s que es el nuevo mínimo para NFW e ISO
 # newi_min_ch2_NFW = [i for i in ch2_inter_NFW_[0] if i[1] == min(list(zip(*ch2_inter_NFW_[0]))[1])][0][0]

 # newi_min_ch2_ISO = [i for i in ch2_inter_ISO_[0] if i[1] == min(list(zip(*ch2_inter_ISO_[0]))[1])][0][0]

 # new_para_min_ch2_NFW = combs_NFW_ra[newi_min_ch2_NFW]
 # new_para_min_ch2_ISO = combs_ISO_ra[newi_min_ch2_ISO]

 # new_para_min_ch2_b = np.concatenate((new_para_min_ch2_NFW,new_para_min_ch2_ISO))

 # print("Los nuevos parámetros para NFW son: {} con un $\chi^2 = {}$"
 #       .format(new_para_min_ch2_NFW, chis_NFW_ra[newi_min_ch2_NFW]))

 # print("Los nuevos parámetros para ISO son: {} con un $\chi^2 = {}$"
 #       .format(new_para_min_ch2_ISO, chis_ISO_ra[newi_min_ch2_ISO]))



 ```

 %% Cell type:code id: tags:

 ``` python
 # new_para_min_ch2_b
 ```

 %% Cell type:code id: tags:

 ``` python
 # print("La probabilidad cumulativa de encontrar un $\chi^2$ menor al nuevo para NFW es: {}".format(scipy.stats.chi2.cdf(chis_NFW_ra[newi_min_ch2_NFW]*(len(df1)-2),len(df1)-2)))
 # print("La probabilidad cumulativa de encontrar un $\chi^2$ menor al nuevo para ISO es: {}".format(scipy.stats.chi2.cdf(chis_ISO_ra[newi_min_ch2_ISO]*(len(df1)-2),len(df1)-2)))
 ```

 %% Cell type:code id: tags:

 ``` python
 # -- gráfica, en functions.py, toca descomentar la línea que guarda la gráfica
 plot_ellipse_confidence_interval(path,para_range_NFW, best_parmeters, ss="ra", i = 0 )
 ```

 %% Output



 %% Cell type:code id: tags:

 ``` python
 # -- gráfica, en functions.py, toca descomentar la línea que guarda la gráfica
 plot_ellipse_confidence_interval(path,para_range_ISO, best_parmeters, ss="ra", i=1 )
 ```

 %% Output



 %% Cell type:code id: tags:

 ``` python
-# -- gráfica, en functions.py, toca descomentar la línea que guarda la gráfica
-#grafica_data_ybestfit(path, df1, best_parmeters, ss="ra", mtype = "new", new_para_min_ch2_b = new_para_min_ch2_b)
-```
-
-%% Cell type:code id: tags:
-
-``` python
 # -- guarda los datos que necesitamos para el overleaf, con suerte es solo copiar y pegar cada linea.
 # -- Toca tambien poner mas bonito la M200, python la da toda fea (como 1.00e+10)

 # file1 = open("results_old.txt","a")
 # file1.write("\n{}&{}&{:3e}&{:2f}&{:2f}&{:2f}&{:2f}&{:2f}&".format(
 # path.split("/")[-1][:-4][5:] # saca el N de la galaxia
 # ,len(df1) # saca la longitud de los datos
 # ,best_parmeters["M200_NFW(Msun)"][0] # aca toca cambiarle el 0 (ra), por 1(r) o 2(a)
 # ,best_parmeters["c_NFW"][0] # aca toca cambiarle el 0 (ra), por 1(r) o 2(a)
-# ,chi2_NFW_bestfit # chi2 para NFW con los parámetros que da leastsq
+# ,chi2_NFW_bestfit_red # chi2 para NFW con los parámetros que da leastsq
 # ,abs(best_parmeters["rho0_ISO(Msun/pc3)"][0]) # aaca toca cambiarle el 0 (ra), por 1(r) o 2(a)
 # ,abs(best_parmeters["Rc_ISO(pc)"][0]) # aca toca cambiarle el 0 (ra), por 1(r) o 2(a)
-# ,chi2_ISO_bestfit)) # chi2 para ISO con los parámetros que da leastsq
+# ,chi2_ISO_bestfit_red)) # chi2 para ISO con los parámetros que da leastsq
 # file1.close()
 ```

 %% Cell type:code id: tags:

 ``` python
 ```

--- a/notebook/confidence_interval_ellipse1.ipynb
+++ b/notebook/confidence_interval_ellipse1.ipynb