No estoy buscando un método plug and play como BEST en R, sino más bien una explicación matemática de algunos de los métodos bayesianos que puedo usar para probar la diferencia entre la media de dos muestras.
Respuestas
¿Demasiados anuncios?La excelente respuesta del usuario 1068430 implementada en Python
import numpy as np
from pylab import plt
def dnorm(x, mu, sig):
return 1/(sig * np.sqrt(2 * np.pi)) * np.exp(-(x - mu)**2 / (2 * sig**2))
def dexp(x, l):
return l * np.exp(- l*x)
def like(parameters):
[mu1, sig1, mu2, sig2] = parameters
return dnorm(sample1, mu1, sig1).prod()*dnorm(sample2, mu2, sig2).prod()
def prior(parameters):
[mu1, sig1, mu2, sig2] = parameters
return dnorm(mu1, pooled.mean(), 1000*pooled.std()) * dnorm(mu2, pooled.mean(), 1000*pooled.std()) * dexp(sig1, 0.1) * dexp(sig2, 0.1)
def posterior(parameters):
[mu1, sig1, mu2, sig2] = parameters
return like([mu1, sig1, mu2, sig2])*prior([mu1, sig1, mu2, sig2])
#create samples
sample1 = np.random.normal(100, 3, 8)
sample2 = np.random.normal(100, 7, 10)
pooled= np.append(sample1, sample2)
plt.figure(0)
plt.hist(sample1)
plt.hold(True)
plt.hist(sample2)
plt.show(block=False)
mu1 = 100
sig1 = 10
mu2 = 100
sig2 = 10
parameters = np.array([mu1, sig1, mu2, sig2])
niter = 10000
results = np.zeros([niter, 4])
results[1,:] = parameters
for iteration in np.arange(2,niter):
candidate = parameters + np.random.normal(0,0.5,4)
ratio = posterior(candidate)/posterior(parameters)
if np.random.uniform() < ratio:
parameters = candidate
results[iteration,:] = parameters
#burn-in
results = results[499:niter-1,:]
mu1 = results[:,1]
mu2 = results[:,3]
d = (mu1 - mu2)
p_value = np.mean(d > 0)
plt.figure(1)
plt.hist(d,normed = 1)
plt.show()
Con un análisis bayesiano, tiene más cosas que especificar (que en realidad es algo bueno, ya que brinda mucha más flexibilidad y capacidad para modelar lo que cree que es la verdad). ¿Está asumiendo normales para las probabilidades? ¿Los 2 grupos tendrán la misma varianza?
Un enfoque sencillo es modelar las 2 medias (y 1 o 2 varianzas / dispersiones) y luego mirar el posterior en la diferencia de las 2 medias y / o el Intervalo creíble en la diferencia de las 2 medias.