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)


Las medias muestrales son
1.99758902456 0.995554765082
Las desviaciones estandar son
0.0100175542603 0.0199419979559

In [158]:
plt.plot(sample[0,:],sample[1,:],'.')


Out[158]:
[<matplotlib.lines.Line2D at 0x28127b079e8>]

:D Exito. Os buscais la vida con ponerlo en python 2


In [ ]: