sistemas - runge kutta 4 orden python
Ajustar una función no lineal a datos/observaciones con pyMCMC/pyMC (2)
Estoy tratando de ajustar algunos datos con una función (es) gaussiana (y más compleja). He creado un pequeño ejemplo a continuación.
Mi primera pregunta es, ¿lo estoy haciendo bien?
Mi segunda pregunta es, ¿cómo agrego un error en la dirección x, es decir, en la posición x de las observaciones / datos?
Es muy difícil encontrar buenas guías sobre cómo hacer este tipo de regresión en pyMC. Tal vez debido a que es más fácil usar algunos mínimos cuadrados, o un enfoque similar, al final tengo muchos parámetros y necesito ver qué tan bien podemos restringirlos y comparar diferentes modelos, pyMC pareció la mejor opción para eso.
import pymc
import numpy as np
import matplotlib.pyplot as plt; plt.ion()
x = np.arange(5,400,10)*1e3
# Parameters for gaussian
amp_true = 0.2
size_true = 1.8
ps_true = 0.1
# Gaussian function
gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps
f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true )
# add noise to the data points
noise = np.random.normal(size=len(x)) * .02
f = f_true + noise
f_error = np.ones_like(f_true)*0.05*f.max()
# define the model/function to be fitted.
def model(x, f):
amp = pymc.Uniform(''amp'', 0.05, 0.4, value= 0.15)
size = pymc.Uniform(''size'', 0.5, 2.5, value= 1.0)
ps = pymc.Normal(''ps'', 0.13, 40, value=0.15)
@pymc.deterministic(plot=False)
def gauss(x=x, amp=amp, size=size, ps=ps):
e = -1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.))
return amp*np.exp(e)+ps
y = pymc.Normal(''y'', mu=gauss, tau=1.0/f_error**2, value=f, observed=True)
return locals()
MDL = pymc.MCMC(model(x,f))
MDL.sample(1e4)
# extract and plot results
y_min = MDL.stats()[''gauss''][''quantiles''][2.5]
y_max = MDL.stats()[''gauss''][''quantiles''][97.5]
y_fit = MDL.stats()[''gauss''][''mean'']
plt.plot(x,f_true,''b'', marker=''None'', ls=''-'', lw=1, label=''True'')
plt.errorbar(x,f,yerr=f_error, color=''r'', marker=''.'', ls=''None'', label=''Observed'')
plt.plot(x,y_fit,''k'', marker=''+'', ls=''None'', ms=5, mew=2, label=''Fit'')
plt.fill_between(x, y_min, y_max, color=''0.5'', alpha=0.5)
plt.legend()
Me doy cuenta de que podría tener que ejecutar más iteraciones, usar quemado y adelgazamiento al final. La figura que representa los datos y el ajuste se ve a continuación.
Las cifras de pymc.Matplot.plot (MDL) se ven así, mostrando distribuciones muy altas. Esto es bueno, ¿verdad?
Mi primera pregunta es, ¿lo estoy haciendo bien?
¡Sí! Es necesario incluir un período de quemado, que usted sabe. Me gusta tirar la primera mitad de mis muestras. No necesita hacer ningún adelgazamiento, pero a veces hará que su trabajo posterior al MCMC sea más rápido de procesar y más pequeño de almacenar.
La única otra cosa que aconsejo es establecer una semilla aleatoria, para que sus resultados sean "reproducibles": np.random.seed(12345) hará el truco.
Ah, y si realmente estuviera dando demasiados consejos, diría que import seaborn para hacer que los resultados de matplotlib un poco más hermosos.
Mi segunda pregunta es, ¿cómo agrego un error en la dirección x, es decir, en la posición x de las observaciones / datos?
Una forma es incluir una variable latente para cada error. Esto funciona en su ejemplo, pero no será factible si tiene muchas más observaciones. Te daré un pequeño ejemplo para que comiences por este camino:
# add noise to observed x values
x_obs = pm.rnormal(mu=x, tau=(1e4)**-2)
# define the model/function to be fitted.
def model(x_obs, f):
amp = pm.Uniform(''amp'', 0.05, 0.4, value= 0.15)
size = pm.Uniform(''size'', 0.5, 2.5, value= 1.0)
ps = pm.Normal(''ps'', 0.13, 40, value=0.15)
x_pred = pm.Normal(''x'', mu=x_obs, tau=(1e4)**-2) # this allows error in x_obs
@pm.deterministic(plot=False)
def gauss(x=x_pred, amp=amp, size=size, ps=ps):
e = -1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.))
return amp*np.exp(e)+ps
y = pm.Normal(''y'', mu=gauss, tau=1.0/f_error**2, value=f, observed=True)
return locals()
MDL = pm.MCMC(model(x_obs, f))
MDL.use_step_method(pm.AdaptiveMetropolis, MDL.x_pred) # use AdaptiveMetropolis to "learn" how to step
MDL.sample(200000, 100000, 10) # run chain longer since there are more dimensions
Parece que puede ser difícil obtener buenas respuestas si tiene ruido en x e y :
Aquí hay un cuaderno que recoge todo esto .
EDIT: nota importante Esto me ha estado molestando desde hace un tiempo. Las respuestas dadas por mí y Abraham aquí son correctas en el sentido de que agregan variabilidad a x. SIN EMBARGO: tenga en cuenta que no puede simplemente agregar incertidumbre de esta manera para cancelar los errores que tiene en sus valores de x, de modo que regrese a la "x verdadera". Los métodos en esta respuesta pueden mostrarle cómo agregar errores a x afecta su regresión si tiene la verdadera x. Si tienes una x no medida, estas respuestas no te ayudarán. Tener errores en los valores de x es un problema muy difícil de resolver, ya que conduce a una "atenuación" y un "efecto de errores en las variables". La versión corta es: tener errores aleatorios e imparciales en x conduce a un sesgo en sus estimaciones de regresión. Si tiene este problema, consulte Carroll, RJ, Ruppert, D., Crainiceanu, CM y Stefanski, LA, 2006. Error de medición en modelos no lineales: una perspectiva moderna . Chapman y Hall / CRC., O para un enfoque bayesiano, Gustafson, P., 2003. Error de medición y error de clasificación en estadísticas y epidemiología: impactos y ajustes bayesianos . Prensa CRC. Terminé resolviendo mi problema específico utilizando el método SIMEX de Carroll et al. Junto con PyMC3. Los detalles se encuentran en Carstens, H., Xia, X. y Yadavalli, S., 2017. Método de calibración de medidor de energía de bajo costo para medición y verificación. Energía aplicada, 188, pp.563-575. También está disponible en ArXiv.
Convertí la respuesta de Abraham Flaxman en PyMC3, en caso de que alguien la necesite. Algunos cambios muy menores, pero pueden ser confusos sin embargo.
La primera es que el decorador determinista @Deterministic se reemplaza por una función de llamada similar a la distribución var=pymc3.Deterministic() . Segundo, cuando se genera un vector de variables aleatorias normalmente distribuidas,
rvs = pymc2.rnormal(mu=mu, tau=tau)
es reemplazado por
rvs = pymc3.Normal(''var_name'', mu=mu, tau=tau,shape=size(var)).random()
El código completo es el siguiente:
import numpy as np
from pymc3 import *
import matplotlib.pyplot as plt
# set random seed for reproducibility
np.random.seed(12345)
x = np.arange(5,400,10)*1e3
# Parameters for gaussian
amp_true = 0.2
size_true = 1.8
ps_true = 0.1
#Gaussian function
gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps
f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true )
# add noise to the data points
noise = np.random.normal(size=len(x)) * .02
f = f_true + noise
f_error = np.ones_like(f_true)*0.05*f.max()
with Model() as model3:
amp = Uniform(''amp'', 0.05, 0.4, testval= 0.15)
size = Uniform(''size'', 0.5, 2.5, testval= 1.0)
ps = Normal(''ps'', 0.13, 40, testval=0.15)
gauss=Deterministic(''gauss'',amp*np.exp(-1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.)))+ps)
y =Normal(''y'', mu=gauss, tau=1.0/f_error**2, observed=f)
start=find_MAP()
step=NUTS()
trace=sample(2000,start=start)
# extract and plot results
y_min = np.percentile(trace.gauss,2.5,axis=0)
y_max = np.percentile(trace.gauss,97.5,axis=0)
y_fit = np.percentile(trace.gauss,50,axis=0)
plt.plot(x,f_true,''b'', marker=''None'', ls=''-'', lw=1, label=''True'')
plt.errorbar(x,f,yerr=f_error, color=''r'', marker=''.'', ls=''None'', label=''Observed'')
plt.plot(x,y_fit,''k'', marker=''+'', ls=''None'', ms=5, mew=1, label=''Fit'')
plt.fill_between(x, y_min, y_max, color=''0.5'', alpha=0.5)
plt.legend()
Lo que resulta en
Para errores en x (note el sufijo ''x'' a las variables):
# define the model/function to be fitted in PyMC3:
with Model() as modelx:
x_obsx = pm3.Normal(''x_obsx'',mu=x, tau=(1e4)**-2, shape=40)
ampx = Uniform(''ampx'', 0.05, 0.4, testval=0.15)
sizex = Uniform(''sizex'', 0.5, 2.5, testval=1.0)
psx = Normal(''psx'', 0.13, 40, testval=0.15)
x_pred = Normal(''x_pred'', mu=x_obsx, tau=(1e4)**-2*np.ones_like(x_obsx),testval=5*np.ones_like(x_obsx),shape=40) # this allows error in x_obs
gauss=Deterministic(''gauss'',ampx*np.exp(-1*(np.pi**2*sizex*x_pred/(3600.*180.))**2/(4.*np.log(2.)))+psx)
y = Normal(''y'', mu=gauss, tau=1.0/f_error**2, observed=f)
start=find_MAP()
step=NUTS()
tracex=sample(20000,start=start)
Lo que resulta en:
La última observación es que al hacer.
traceplot(tracex[100:])
plt.tight_layout();
(resultado no mostrado), podemos ver que sizex parece estar sufriendo de "atenuación" o "dilución de regresión" debido al error en la medición de x .