Gibbs-Sampler:
Sea $D=(x_i,y_i)_{i=1,..,N}$ una muestra de una distribucion gaussiana 2D. Sean $\mu_x$, $\mu_y$, $\sigma_x$, $\sigma_y$ y $\rho$ las correspondientes medias, desviaciones estandar y el coeficiente de correlacion.
Por el ejercicio 3.II sabemos que el posterior es:
$$P(\mu_x,\mu_y|D)=\frac{c}{(2\pi\sigma_x\sigma_y\sqrt{1-\rho^2})^N}e^{-\frac{1}{2}\sum_{i=1}^N{A_i}}$$Donde $A_i$ viene dada por:
$$A_i=\frac{(x_i-\mu_x)^2}{\sigma_x^2(1-\rho^2)}+\frac{(y_i-\mu_y)^2}{\sigma_y^2(1-\rho^2)}-\frac{2\rho}{\sigma_x\sigma_y(1-\rho^2)}(y_i-\mu_y)(x_i-\mu_x)$$Desarrollando $A_i$ se llega a que:
$$A_i=\frac{x^2_i+\mu_x^2-2\mu_x x_i}{\sigma_x^2(1-\rho^2)}+\frac{(y_i-\mu_y)^2}{\sigma_y^2(1-\rho^2)}-\frac{2\rho}{\sigma_x\sigma_y(1-\rho^2)}x_i y_i+\frac{2\rho}{\sigma_x\sigma_y(1-\rho^2)}y_i\mu_x+\frac{2\rho}{\sigma_x\sigma_y(1-\rho^2)}x_i\mu_y-\frac{2\rho}{\sigma_x\sigma_y(1-\rho^2)}\mu_x\mu_y$$Por lo tanto se tiene que:
$$P(\mu_x | \mu_y;D)\propto \exp\bigg(-\frac{1}{2}\sum_{i=1}^N\Big(\frac{\mu_x^2-2\mu_x x_i}{\sigma_x^2(1-\rho^2)}+\frac{2\rho}{\sigma_x\sigma_y(1-\rho^2)}y_i\mu_x-\frac{2\rho}{\sigma_x\sigma_y(1-\rho^2)}\mu_x\mu_y\Big)\bigg)$$Esto es debido a que $P(\mu_x | \mu_y;D)$ es la probabilidad de obtener $\mu_x$ cuando $\mu_y$ esta FIJO, por lo tanto cualquier termino donde no haya un $\mu_x$ sera constante.
Desarrollando aun mas esta expresion
$$P(\mu_x | \mu_y;D)\propto \exp\bigg(-\frac{1}{2}\sum_{i=1}^N{\frac{\mu_x^2\sigma_y-2\mu_x x_i\sigma_y+2\rho\sigma_x\mu_x(y_i-\mu_y)}{\sigma_x^2\sigma_y(1-\rho^2)}}\bigg)=\exp\bigg(-\frac{1}{2}\sum_{i=1}^N{\frac{\mu_x^2-\mu_x\big(2x_i-2\rho\frac{\sigma_x}{\sigma_y}(y_i-\mu_y)\big)}{\sigma_x^2(1-\rho^2)}}\bigg)$$Entonces, ya que $\mu_j$ (j=x,y) no depende del sumatorio se tiene que:
$$\sum_{i=1}^N{\Big(\mu_x^2-2\mu_x\big(x_i-\rho\frac{\sigma_x}{\sigma_x}(y_i-\mu_y)\big)\Big)}=\mu_x^2 N-2\mu_x\bigg(\sum_{i=1}^N{x_i}-\rho\frac{\sigma_x}{\sigma_y}\Big(\sum_{i=1}^N{y_i}-N\mu_y\Big)\bigg)=N\Big(\mu_x^2-2\mu_x\big(\frac{1}{N}\sum_{i=1}^N{x_i}-\rho\frac{\sigma_x}{\sigma_y}(\frac{1}{N}\sum_{i=1}^N{y_i}-\mu_y)\big)\Big)$$Definiendo las medias muestrales como $$\bar{x}=\frac{1}{N}\sum_{i=1}^N{x}_i \quad ; \quad \bar{y}=\frac{1}{N}\sum_{i=1}^N{y}_i$$ Se obtiene la expresion: $$N\Big(\mu_x^2-2\mu_x\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)\Big)$$ Y por lo tanto: $$P(\mu_x | \mu_y;D)\propto \exp\bigg(-\frac{N}{2\sigma_x^2(1-\rho^2)}\Big(\mu_x^2-2\mu_x\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)\Big)\bigg)$$
Sumando y restando el termino $\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)^2$:
$$P(\mu_x | \mu_y;D)\propto \exp\bigg(-\frac{N}{2\sigma_x^2(1-\rho^2)}\Big(\mu_x^2-2\mu_x\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)+\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)^2-\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)^2\Big)\bigg)$$Y haciendo cuadrados: $$P(\mu_x | \mu_y;D)\propto \exp\Bigg(-\frac{N}{2\sigma_x^2(1-\rho^2)}\bigg(\Big(\mu_x-\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)\Big)^2-\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)^2\bigg)\Bigg) \propto \exp\bigg(-\frac{\Big(\mu_x-\big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)\big)\Big)^2}{2\frac{\sigma_x^2(1-\rho^2)}{N}}\bigg)$$
Que se corresponde con la forma de una distribucion normal con media $\mu=\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y)$ y desviacion estandar $\sigma=\frac{\sigma_x\sqrt{1-\rho^2}}{\sqrt{N}}$.
Se concluye que: $$P(\mu_x | \mu_y;D)\equiv \mathcal{N}\Big(\bar{x}-\rho\frac{\sigma_x}{\sigma_y}(\bar{y}-\mu_y),\frac{\sigma_x\sqrt{1-\rho^2}}{\sqrt{N}}\Big)$$ $$P(\mu_y | \mu_x;D)\equiv \mathcal{N}\Big(\bar{y}-\rho\frac{\sigma_y}{\sigma_x}(\bar{x}-\mu_x),\frac{\sigma_y\sqrt{1-\rho^2}}{\sqrt{N}}\Big)$$ Donde la segunda probabilidad condicionada se obtiene directamente debido a la simetria que existe entre $\mu_x$ y $\mu_y$ en la ecuacion original.
In [2]:
%matplotlib inline
import numpy as np
import math
import scipy.stats as stat # extra statistical functions (the basic are included in numpy)
import scipy.optimize as opt # optimization and root finding package
import scipy.misc as msc
import scipy.special as spc
import matplotlib.pyplot as plt
In [3]:
N_3I = 10**4
mu1 = 2
mu2 = 1
sigma1 = 1
sigma2 = 2
rho = 0.80
y1 = np.random.normal(size=N_3I)
y2 = np.random.normal(size=N_3I)
x1 = mu1 + sigma1*y1
x2 = mu2 + rho*sigma2*y1 + np.sqrt(1-rho**2)*sigma2*y2
D = np.transpose(np.array([x1,x2]))
Mean1=np.mean(D[:,0])
Mean2=np.mean(D[:,1])
In [157]:
N = 10**4
Sig1=sigma1*np.sqrt(1-rho**2)/np.sqrt(len(D))
Sig2=sigma2*np.sqrt(1-rho**2)/np.sqrt(len(D))
sample = np.zeros((2,N+1))
D1=Mean1-rho*sigma1/sigma2*Mean2
D2=Mean2-rho*sigma2/sigma1*Mean1
for i in range(N):
sample[0,i+1] = np.random.normal(D1+rho*sigma1/sigma2*sample[1,i], Sig1)
sample[1,i+1] = np.random.normal(D2+rho*sigma2/sigma1*sample[0,i+1], Sig2)
SampleMean1=np.mean(sample[0,50:])
SampleMean2=np.mean(sample[1,50:])
print('Las medias muestrales son')
print(SampleMean1, SampleMean2)
SampleStd1=np.std(sample[0,50:])
SampleStd2=np.std(sample[1,50:])
print('Las desviaciones estandar son')
print(SampleStd1, SampleStd2)
In [158]:
plt.plot(sample[0,:],sample[1,:],'.')
Out[158]:
:D Exito. Os buscais la vida con ponerlo en python 2
In [ ]: