Pruebas al sistema de creación de distribuciones PyMC.

Matemáticas/Distribuciones.py

@@ -106,6 +106,7 @@
                'límites': (0, np.inf),
                'tipo': 'cont'
                },
+
          'BirnbaumSaunders': {'scipy': estad.fatiguelife,
                               'pymc': None,
                               'límites': (0, np.inf),
@@ -329,11 +330,15 @@
                         'límites': (0, np.inf),
                         'tipo': 'cont'
                         },
-         'TNoCentral': {'scipy': estad.nct,
-                        'pymc': NoncentralT,
-                        'límites': (-np.inf, np.inf),
-                        'tipo': 'cont'
-                        },
+
+         # Desactivada por complicación de conversión PyMC-SciPy
+         #
+         # 'TNoCentral': {'scipy': estad.nct,
+         #                'pymc': NoncentralT,
+         #                'límites': (-np.inf, np.inf),
+         #                'tipo': 'cont'
+         #                },
+
          'Normal': {'scipy': estad.norm,
                     'pymc': Normal,
                     'límites': (-np.inf, np.inf),
@@ -395,11 +400,14 @@
                           'límites': (-1, 1),
                           'tipo': 'cont'
                           },
-         'T': {'scipy': estad.t,
-               'pymc': T,
-               'límites': (-np.inf, np.inf),
-               'tipo': 'cont'
-               },
+
+         # Desactivada por complicación de conversión PyMC-SciPy
+         # 'T': {'scipy': estad.t,
+         #       'pymc': T,
+         #       'límites': (-np.inf, np.inf),
+         #       'tipo': 'cont'
+         #       },
+
          'Triang': {'scipy': estad.triang,
                     'pymc': None,
                     'límites': (0, 1),  # El límite es (a, b)
@@ -411,11 +419,15 @@
                               'límites': (0, 1),  # El límite es (0, b)
                               'tipo': 'cont'
                               },
-         'NormalTrunc': {'scipy': estad.truncnorm,
-                         'pymc': TruncatedNormal,
-                         'límites': (0, 1),  # El límite es (a, b)
-                         'tipo': 'cont'
-                         },
+
+         # Desactivada por falta de posibilidad de fijar sus límites para la optimización.
+         #
+         # 'NormalTrunc': {'scipy': estad.truncnorm,
+         #                 'pymc': TruncatedNormal,
+         #                 'límites': (0, 1),
+         #                 'tipo': 'cont'
+         #                 },
+
          'TukeyLambda': {'scipy': estad.tukeylambda,
                          'pymc': None,
                          'límites': (-np.inf, np.inf),
@@ -426,11 +438,14 @@
                       'límites': (0, 1),  # El límite es (a, b)
                       'tipo': 'cont'
                       },
-         'VonMises': {'scipy': estad.vonmises,
-                      'pymc': VonMises,
-                      'límites': (-mat.pi, mat.pi),
-                      'tipo': 'cont'
-                      },
+
+         # Desactivada por complicación de conversión PyMC-SciPy
+         # 'VonMises': {'scipy': estad.vonmises,
+         #              'pymc': VonMises,
+         #              'límites': (-mat.pi, mat.pi),
+         #              'tipo': 'cont'
+         #              },
+
          'VonMisesLín': {'scipy': estad.vonmises_line,
                          'pymc': None,
                          'límites': (-mat.pi, mat.pi),

Matemáticas/Ecuaciones.py

@@ -25,7 +25,7 @@
                                                                    'inter': 'presa'}},
                                          'Logístico Depredación': {'K': {'límites': (0, np.inf),
                                                                    'inter': 'presa'}},
-                                         'Externa cultivo': {}
+                                         'Externo Cultivo': {}
                                          }
                             },
 
@@ -96,9 +96,7 @@
                                          'Kovai': {'a': {'límites': (0, np.inf),
                                                          'inter': 'presa'},
                                                    'b': {'límites': (0, 1),
-                                                         'inter': None},
-                                                   'c': {'límites': (1, np.inf),
-                                                         'inter': None}
+                                                         'inter': 'presa'},
                                                    }
                                          },
                             },

Matemáticas/NuevoIncert.py

@@ -216,17 +216,14 @@ def dist_a_texto(dist):
     return texto_dist
 
 
-def ajustar_dist(datos, límites, cont, usar_pymc=False, nombre=None):
+def ajustar_dist(datos, límites, cont, usar_pymc=False, nombre=None, lista_dist=None):
     """
     Esta función, tomando las límites teoréticas de una distribución y una serie de datos proveniendo de dicha
       distribución, escoge la distribución de Scipy o PyMC la más apropriada y ajusta sus parámetros.
 
     :param datos: Un vector de valores del parámetro
     :type datos: np.ndarray
 
-    :param nombre: El nombre del variable, si vamos a generar un variable de PyMC
-    :type nombre: str
-
     :param cont: Determina si la distribución es contínua (en vez de discreta)
     :type cont: bool
 
@@ -236,6 +233,12 @@ def ajustar_dist(datos, límites, cont, usar_pymc=False, nombre=None):
     :param límites: Las límites teoréticas de la distribucion (p. ej., (0, np.inf), (-np.inf, np.inf), etc.)
     :type límites: tuple
 
+    :param nombre: El nombre del variable, si vamos a generar un variable de PyMC
+    :type nombre: str
+
+    :param lista_dist: Una lista de los nombres de distribuciones a considerar. dist=None las considera todas.
+    :type lista_dist: list
+
     :return: Distribución PyMC o de Scipyy su ajuste (p)
     :rtype: (pymc.Stochastic | estad.rv_frozen, float)
 
@@ -259,6 +262,12 @@ def ajustar_dist(datos, límites, cont, usar_pymc=False, nombre=None):
 
     dists_potenciales = [x for x in Ds.dists if Ds.dists[x]['tipo'] == categ_dist]
 
+    if lista_dist is not None:
+        dists_potenciales = [x for x in dists_potenciales if x in lista_dist]
+
+        if len(dists_potenciales) == 0:
+            raise ValueError('Ninguna de las distribuciones especificadas es apropiada para el tipo de distribución.')
+
     # Si queremos generar una distribución PyMC, guardar únicamente las distribuciones con objeto de PyMC disponible
     if usar_pymc is True:
         dists_potenciales = [x for x in dists_potenciales if Ds.dists[x]['pymc'] is not None]
@@ -311,29 +320,33 @@ def ajustar_dist(datos, límites, cont, usar_pymc=False, nombre=None):
                     restric = {'floc': mín_parám, 'fscale': máx_parám - mín_parám}
 
             # Ajustar los parámetros de la distribución SciPy para caber con los datos.
-            args = dic_dist['scipy'].fit(datos, **restric)
+            if nombre_dist != 'Uniforme':
+                args = dic_dist['scipy'].fit(datos, **restric)
+            else:
+                # Para distribuciones uniformes, no hay nada que calibrar.
+                args = tuple(restric.values())
 
             # Medir el ajuste de la distribución
             p = estad.kstest(rvs=datos, cdf=nombre_scipy, args=args)[1]
 
             # Si el ajuste es mejor que el mejor ajuste anterior...
-            if p > mejor_ajuste['p']:
+            if p >= mejor_ajuste['p']:
 
                 # Guardarlo
                 mejor_ajuste['p'] = p
 
                 # Guardar también el objeto de la distribución, o de PyMC, o de SciPy, según lo que queremos
                 if usar_pymc:
                     # Convertir los argumentos a formato PyMC
-                    args, transform = paráms_scipy_a_pymc(tipo_dist=nombre_dist, paráms=args)
+                    args_pymc, transform = paráms_scipy_a_pymc(tipo_dist=nombre_dist, paráms=args)
 
                     # Y crear la distribución.
-                    dist = dic_dist['pymc'](nombre, *args)
+                    dist = dic_dist['pymc'](nombre, *args_pymc)
 
                     if transform['sum'] != 0:
                         dist = dist + transform['sum']
                     if transform['mult'] != 1:
-                        dist = dist + transform['mult']
+                        dist = dist * transform['mult']
                     mejor_ajuste['dist'] = dist
 
                 else:
@@ -438,6 +451,10 @@ def rango_a_texto_dist(rango, certidumbre, líms, cont):
     mín = líms[0]
     máx = líms[1]
 
+    # Si el rango está al revés, arreglarlo.
+    if rango[0] > rango[1]:
+        rango = (rango[1], rango[0])
+
     # Asegurarse de que el rango cabe en los límites
     if rango[0] < mín or rango[1] > máx:
         raise ValueError('El rango tiene que caber entre los límites teoréticos del variable.')
@@ -450,7 +467,14 @@ def rango_a_texto_dist(rango, certidumbre, líms, cont):
             dist = 'UnifDiscr~({}, {})'.format(rango[0], rango[1])
 
     else:
-        # Si hay incertidumbre, asignar una distribución según cada caso posible
+        # Si hay incertidumbre...
+
+        # Primero, hay que asegurarse que el máximo y mínimo del rango no sean iguales
+        if rango[0] == rango[1]:
+            raise ValueError('Un rango para una distribucion con certidumbre inferior a 1 no puede tener valores'
+                             'mínimos y máximos iguales')
+
+        # Ahora, asignar una distribución según cada caso posible
         if mín == -np.inf:
             if máx == np.inf:  # El caso (-np.inf, np.inf)
                 if cont:
@@ -761,15 +785,15 @@ def paráms_scipy_a_pymc(tipo_dist, paráms):
         transform_pymc['sum'] = paráms[1]
 
     elif tipo_dist == 'MitadCauchy':
-        paráms_pymc = (paráms[0], 1 / paráms[1])
+        paráms_pymc = (paráms[0], paráms[1])
 
     elif tipo_dist == 'MitadNormal':
-        paráms_pymc = (1 / paráms[1], )
+        paráms_pymc = (1 / paráms[1]**2, )
         transform_pymc['sum'] = paráms[0]
 
     elif tipo_dist == 'GammaInversa':
-        paráms_pymc = (paráms[0], paráms[1])
-        transform_pymc['sum'] = paráms[2]
+        paráms_pymc = (paráms[0], paráms[2])
+        transform_pymc['sum'] = paráms[1]
 
     elif tipo_dist == 'Laplace':
         paráms_pymc = (paráms[0], 1 / paráms[1])
@@ -778,7 +802,8 @@ def paráms_scipy_a_pymc(tipo_dist, paráms):
         paráms_pymc = (paráms[0], 1 / paráms[1])
 
     elif tipo_dist == 'LogNormal':
-        paráms_pymc = (0, 1 / (paráms[0]**2))
+        paráms_pymc = (np.log(paráms[2]), 1 / (paráms[0]**2))
+        transform_pymc['mult'] = paráms[2]
         transform_pymc['sum'] = paráms[1]
 
     elif tipo_dist == 'TNoCentral':
@@ -794,10 +819,12 @@ def paráms_scipy_a_pymc(tipo_dist, paráms):
     elif tipo_dist == 'T':
         paráms_pymc = (paráms[0], )
         transform_pymc['sum'] = paráms[1]
-        transform_pymc['mult'] = paráms[2]
+        transform_pymc['mult'] = 1 / np.sqrt(paráms[2])
 
     elif tipo_dist == 'NormalTrunc':
-        paráms_pymc = (paráms[2], 1 / paráms[3]**2, paráms[0], paráms[1])
+        mu, sigma = paráms[2], paráms[3]
+        mín, máx = min(paráms[0], paráms[1]), max(paráms[0], paráms[1])  # SciPy, aparamente, los puede inversar
+        paráms_pymc = (mu, 1 / sigma**2, mín * sigma + mu, máx * sigma + mu)
 
     elif tipo_dist == 'Uniforme':
         paráms_pymc = (paráms[0], paráms[1] + paráms[0])

Matemáticas/Pruebas.py

@@ -0,0 +1,100 @@
+import numpy as np
+import scipy.stats as estad
+import pymc
+import matplotlib.pyplot as dib
+from warnings import warn as avisar
+
+import Matemáticas.NuevoIncert as Inc
+import Matemáticas.Distribuciones as Ds
+
+"""
+Un módulo de pruebas, un poco inelegante.
+"""
+
+# Pruebas para el convertidor de distribuciones SciPy a PyMC en NuevoIncert
+
+n_rep = 10000
+
+dibujar = False
+
+con_errores = []
+
+for nombre, dist in sorted(Ds.dists.items()):
+
+    if dist['tipo'] != 'cont':
+        continue
+
+    if dist['scipy'] is None or dist['pymc'] is None:
+        continue
+
+    mín, máx = dist['límites']
+    if mín == -np.inf:
+        if máx == np.inf:
+            núms = estad.norm.rvs(10, 20, size=100)
+        else:
+            avisar('No se pudo comprobar la distribución %s' % nombre)
+            continue
+
+    else:
+        if máx == np.inf:
+            núms = estad.gamma.rvs(1, loc=mín, size=100)
+        else:
+            núms = estad.uniform.rvs(mín, máx-mín, size=100)
+
+    dist_scipy = Inc.ajustar_dist(datos=núms, límites=dist['límites'], cont=True,
+                                  lista_dist=[nombre])[0]
+    puntos_scipy = dist_scipy.rvs(size=n_rep)
+
+    dist_pymc = Inc.ajustar_dist(datos=núms, límites=dist['límites'], cont=True, usar_pymc=True,
+                                 nombre=nombre, lista_dist=[nombre])[0]
+    if isinstance(dist_pymc, pymc.Stochastic):
+        puntos_pymc = np.array([dist_pymc.rand() for x in range(n_rep)])
+    elif isinstance(dist_pymc, pymc.Deterministic):
+        puntos_pymc = np.empty(n_rep)
+        dist_pariente = min(dist_pymc.extended_parents)
+        for i in range(n_rep):
+            dist_pariente.rand()
+            puntos_pymc[i] = dist_pymc.value
+    else:
+        raise ValueError
+
+    scipy_de_pymc = Inc.ajustar_dist(datos=puntos_pymc, límites=dist['límites'], cont=True,
+                                     lista_dist=[nombre])[0]
+
+    if dibujar:
+
+        dib.hist(núms, normed=True, color='red', histtype='stepfilled', alpha=0.2, bins=100)
+
+        dib.hist(puntos_scipy, normed=True, color='blue', histtype='stepfilled', alpha=0.2, bins=100)
+
+        dib.hist(puntos_pymc, normed=True, color='green', histtype='stepfilled', alpha=0.2, bins=100)
+
+        x = np.linspace(dist_scipy.ppf(0.01), dist_scipy.ppf(0.99), 100)
+        dib.plot(x, dist_scipy.pdf(x), 'b-', lw=2, alpha=0.6)
+        dib.plot(x, scipy_de_pymc.pdf(x), 'g-', lw=2, alpha=0.6)
+
+        dib.title(nombre)
+        dib.show()
+
+    error = 0
+    for n, i in enumerate(scipy_de_pymc.args):
+        if i == 0:
+            e = i - dist_scipy.args[n]
+        else:
+            e = (i - dist_scipy.args[n]) / dist_scipy.args[n]
+        error = max(error, abs(e))
+
+    if error >= 0.05:
+        mensaje = '¡¡¡!!!'
+        con_errores.append(nombre)
+    else:
+        mensaje = 'Todo bien.'
+
+    print(nombre, error, mensaje, dist_scipy.args, scipy_de_pymc.args)
+
+if len(con_errores):
+    print('Verificar:', con_errores)
+else:
+    print('¡Todas las distribuciones pasaron!')
+
+# Arreglar: Errores con T, TNo, etc...