Python: Ajuste MCMC a datos con incertidumbres variables (usando lmfit)

Answer

1. 

   davy.ai

   August 6, 2023 at 10:43 am

   Para ajustar ambos conjuntos de datos simultáneamente teniendo en cuenta las incertidumbres de los datos variables, se puede utilizar un enfoque bayesiano. Esto implica definir una función de verosimilitud que incorpore las incertidumbres, y utilizar métodos MCMC (cadenas de Markov Monte Carlo) para muestrear de la distribución posterior de los parámetros del modelo.

   Una forma de hacer esto en Python es utilizar el paquete emcee en lugar de lmfit. emcee es un paquete especializado en MCMC que se puede utilizar para ajustar modelos complejos a datos.

   Aquí tienes un ejemplo de cómo implementar este enfoque:

   import numpy as np
   import emcee
   
   # Definir el modelo
   def modelo(x, s, d, c):
       return (s/np.pi) * np.arctan(x/d) + c
   
   # Definir la función de log-verosimilitud, teniendo en cuenta las incertidumbres en y1
   def ln_verosimilitud(p, x, y1, y1_err, y2):
       s, d, c = p
       modelo_y1 = modelo(x, s, d, c)
       residuos_y1 = y1 - modelo_y1
       ln_verosimilitud_y1 = -0.5 * np.sum((residuos_y1 / y1_err)**2)
       ln_verosimilitud_y2 = -0.5 * np.sum((y2 - modelo_y1)**2)
       return ln_verosimilitud_y1 + ln_verosimilitud_y2
   
   # Establecer las condiciones previas sobre los parámetros
   def ln_previo(p):
       s, d, c = p
       if -10 < s < 10 and 0 < d < 30 and -10 < c < 10:
           return 0.0
       else:
           return -np.inf
   
   # Definir la función de log-probabilidad posterior
   def ln_posterior(p, x, y1, y1_err, y2):
       ln_previo_val = ln_previo(p)
       if not np.isfinite(ln_previo_val):
           return -np.inf
       else:
           ln_verosimilitud_val = ln_verosimilitud(p, x, y1, y1_err, y2)
           return ln_previo_val + ln_verosimilitud_val
   
   # Configurar el simulador
   nwalkers = 50
   ndim = 3
   p0 = np.random.rand(nwalkers, ndim)
   simulador = emcee.EnsembleSampler(nwalkers, ndim, ln_posterior, args=(x, y1, y1_err, y2))
   
   # Quemar el simulador
   posicion, probabilidad, estado = simulador.run_mcmc(p0, 1000)
   
   # Reiniciar el simulador y correr por más tiempo
   simulador.reset()
   simulador.run_mcmc(posicion, 10000)
   
   # Obtener las muestras posteriores
   muestras = simulador.chain[:, 1000:, :].reshape((-1, ndim))
   
   # Calcular los valores de mejor ajuste e incertidumbres
   s_mejor, d_mejor, c_mejor = np.median(muestras, axis=0)
   s_std, d_std, c_std = np.std(muestras, axis=0)
   
   

   En este ejemplo, x, y1 y y2 son la variable independiente y dos conjuntos de datos, respectivamente. y1_err es el array de incertidumbres para el primer conjunto de datos. La función ln_verosimilitud calcula la log-verosimilitud para ambos conjuntos de datos, ponderando los residuos del primer conjunto de datos por las incertidumbres. La función ln_previo define la distribución de probabilidad previa para los parámetros, y la función ln_posterior combina la verosimilitud y el previo para calcular la log-probabilidad posterior. Luego, se utiliza el objeto emcee.EnsembleSampler para muestrear de la distribución posterior, utilizando un período de quemado para descartar cualquier muestra inicial deficiente.

   Una vez finalizado el muestreo, se extraen las muestras posteriores y se utilizan para calcular los valores de mejor ajuste e incertidumbres.

Comments are closed.