Hello, nbpresent!

In [ ]:
import nbpresent
nbpresent.__version__

Episodio IV

Willy Pregliasco y Mariano Gómez Berisso, Ago 2016

Distribuciones

Vamos a ver tres distribuciones que aparecen mucho en la Física y en el análisis de los experimentos.

Las tres son distribuciones de una sola variable.

En todos los casos, las distribuciones dependen de unos parámetros $\theta$ y tenemos una serie de $N$ datos medidos $D=\{x_i\}$.

Nos va a interesar la estimación de los parámetros en base a los datos medidos. Lo podemos pensar como un problema estándar de fiteo. Como siempre, aplicamos el teorema de Bayes:

$$ P(\theta \ | D I) = \frac{P(D\ | \theta I)}{P(D\ | I)} \ P(\theta\ | I) $$

En este apunte, vamos a suponer que la prior $P(\theta\ | I)$ es uniforme. El denominador (como siempre) es un factor de normalización que no depende de los parámetros, y por último, las mediciones son independientes. Las dos hipótesis nos permiten escribir:

$$ P(\theta \ | D I) \propto \prod_i P(x_i\ | \theta I) $$

1. Distribución Normal (de Gauss)

La distribución gaussiana es la que más aparece en los experimentos y en la imaginación de los científicos. Hay carreas enteras edificadas en esta única distribución.

Es una distribución simétrica que tiene dos parámetros: un centro $\mu$ y un factor de escala $\sigma$.

La densidad de probabilidad tiene la forma:

$$ G(x,\mu,\sigma) = \frac{1}{\sqrt{2\pi \sigma^2}} \ e^{-\frac{1}{2\sigma^2} (x-\mu)^2} $$

El valor medio coincide con el máximo y es $\mu$. La varianza de la distribución es $\sigma^2$.

Tiene un lugar especial dentro de las distribuciones, porque existe un resultado muy general:

Teorema Central del Límite
Si tengo muchas variables aleatorias que tienen definido un promedio y una varianza, la suma de esas variables es también una variable aleatoria que tiene una distribución que tiende a la gaussiana.

Las consecuencias de este teorema son:

  • si en una medición hay muchos factores que intervienen, la influencia de todos juntos tiende a parecerse a una gaussiana.
  • la combinación de variables gaussianas, también es gaussiana.
  • puedo generar una distribución gaussiana sumando variables con distribuciones muy diferentes.

1.1 Dibujo

In [1]:
### Dibujito de la Gaussiana
%pylab inline
from scipy.stats import norm

figure(figsize=(10,7))
x=linspace(-20,20,200)
title('Distribución Gaussiana con $\mu=5$',fontsize=20)
xlabel('$x$',fontsize=20)
ylabel('$G(x,\mu,\sigma)$',fontsize=20)
plot(x,norm.pdf(x,loc=5,scale=1),label=r'$\sigma=1$',lw=2)
plot(x,norm.pdf(x,loc=5,scale=2),label=r'$\sigma=2$',lw=2)
plot(x,norm.pdf(x,loc=5,scale=5),label=r'$\sigma=5$',lw=2)
grid()
legend();

Populating the interactive namespace from numpy and matplotlib

1.2 Percentiles

Como la gaussiana se utiliza mucho, hay algunas integrales de probabilidad que conviene aprenderse de memoria.

In [2]:
figure(figsize=(10,7))
title('Percentiles en la Gaussiana', fontsize=20)

x=linspace(-5,5,200)
plot(x,norm.pdf(x),lw=2)

x1=linspace(-1,1,100)
fill_between(x1,0,norm.pdf(x1),alpha=.8,color='b')


x2=linspace(-2,2,100)
fill_between(x2,0,norm.pdf(x2),alpha=.4,color='b')

x3=linspace(-3,3,100)
fill_between(x3,0,norm.pdf(x3),alpha=.3,color='k')

text(0  ,0.15,'68%',fontsize=20,ha='center',color='w')
text(1.5,0.03,'95%',fontsize=20,ha='center',color='w')
text(2.9,0.03,'99.7%',fontsize=20,ha='center',color='k')

ylabel('Pdf',fontsize=20)
#xlabel('$x$',fontsize=20)
ylim(0,.45)
xlim(-5,5)
xticks([-3,-2,-1,0,1,2,3],
       ['$-3\sigma$','$-2\sigma$','$\sigma$', '$\mu$', '$\sigma$', '$2\sigma$', '$3\sigma$'], fontsize=20)
grid()

1.3 Como límite de lanzar una moneda

Vamos a usar una variable aleatoria que resulta de lanzar una moneda. Vale 1 si sale cara y 0 si sale ceca. Ambas tienen la misma probabilidad.

Generamos $N$ datos y los sumamos. Hacemos lo mismo $M$ veces. Dibujamos el histograma resultante y comparamos con la gaussiana.

In [3]:
def generaG(N,M):
    '''Sumo N tiradas de una moneda
       Calculo esa variable M veces
    '''
    out=[]
    ## Veces que calculo la variable suma
    for i in range(M):
        ## Tiro N veces la moneda
        x = rand(N)
        x[x>=0.5]=1
        x[x< 0.5]=0
        ## Sumo
        out.append(sum(x))
    return array(out)

def plotH(N,M,**opt):
        '''Dibuja el resultado
        '''
        data = generaG(N,M)
        hist(data, 
             bins=opt['bins'], 
             normed=True,
             histtype='stepfilled',
             alpha=.4,
             color=opt['color']    )

        mu,sigma = mean(data), std(data)
        x = linspace(0,35,200)
        plot(x,norm.pdf(x,loc=mu, scale=sigma),
             color=opt['color'], 
             lw=2,
             label='$N$ = {:d}'.format(N))

figure(figsize=(15,6))        
title('Suma de $N$ tiradas de moneda, comparados con Gaussiana', fontsize=20) 
xlabel('$\sum_i^N \ x_i$', fontsize=15)
ylabel('frecuencia', fontsize=15)
plotH( 3,10000,bins= 4,color='b')        
plotH(15,10000,bins= 7,color='r')
plotH(40,10000,bins=12,color='g')

ylim(0,.6)
xlim(0,30)
legend(fontsize=20)
grid()

1.3 Estimación de Parámetros

Dado un conjunto de datos, tenemos que estimar dos parámetros. Como vimos en la introducción, con priors uniformes resulta que:

$$ P(\mu \sigma \ | D I) \propto \prod_i G(x_i ,\mu, \sigma) = \frac{1}{(\sqrt{2\pi \sigma^2})^N} \ \exp -\frac{1}{2\sigma^2}\ \sum_i (x_i-\mu)^2 $$

Podemos hacer algo de cosmética con la sumatoria:

\begin{eqnarray} \sum_i (x_i-\mu)^2 & = & \sum_i (\mu^2 - 2x_i\mu + x_i^2 ) \\ & = & N \left(\mu^2 -2\mu \overline{x} + \overline{x^2} \right) \\ & = & N \left(\mu^2 -2\mu \overline{x} + \overline{x}^2 - \overline{x}^2 + \overline{x^2} \right) \\ & = & N \left[\ (\mu-\overline{x})^2 + \mathrm{Var}(x)\ \right] \end{eqnarray}

y así llegamos al bonito resultado:

$$ P(\mu \sigma \ | D I) \propto \sigma^{-N} \ \exp -\frac{N}{2\sigma^2}\ \left[\ (\mu-\overline{x})^2 + \mathrm{Var}(x)\ \right] $$
Sacar conclusiones de esta expresión es un trabajo que está lleno de detalles engorrosos. La dificultad no es insalvable, pero resulta poco didáctica en este momento. El problema es que para obtener la distribución de $\mu$ tenemos que marginalizar en $\sigma$, haciendo: $$ P(\mu \ | D I) = \int P(\mu \sigma \ | D I) P(\sigma \ | I) d\sigma $$ y aquí son importantes las priors, que hay que considerarlas con más detalle. No hay grandes profundidades conceptuales, así que vamos a optar por calcular la distribución de $\mu$ suponiendo que conozco $\sigma$ y viceversa. Eso va a estar tan cerca de resultado correcto que va a ser suficiente por ahora.

Distribución de $\mu$ dado $\sigma$

Esta distribución, como función de $\mu$, es una gaussiana. El factor $N$ que aparece en el exponente, hace que la dispersión del estimador de $\mu$ sea más chico en un factor (famoso) de $1/\sqrt{N}$:

\begin{eqnarray} \hat{\mu} &=& \overline{x} \\ \sigma_{\mu} &=& \frac{\sigma}{\sqrt{N}} \end{eqnarray}

Distribución de $\sigma$ dado $\mu$

La distribución como función de $\sigma$ ya no se parece a nada conocido. Por ahora sólo nos interesa el máximo y algunas propiedades generales, de manera que lo vamos a analizar sólo en la línea definida por $\mu=\overline{x}$.

$$ P(\sigma \ |\ \mu=\overline{x}, \ D I) \propto \sigma^{-N} \ \exp -\frac{N}{2\sigma^2}\ \mathrm{Var}(x) $$

La dependencia se encuentra en el factor multiplicativo $\sigma^{-N}$ y el denominador del exponente. Ambos hacen nula la función en $+\infty$ y $0$, respectivamente, así que debe haber un máximo en algún lado.
Como se usa mucho, tiene un nombre, pero como es una función fea, tiene un nombre feo: es la distribución $\chi^2$ inversa. Es asimétrica y se vuelve más concentrada al aumentar $N$.

Si buscamos el máximo (es fácil, hay que hacer cero la derivada) y calculamos la varianza (sale de la Wikipedia) resulta que:

\begin{eqnarray} \hat{\sigma} &=& \sqrt{Var(x)} \\ \sigma_{\sigma} &=& \frac{1}{\sqrt{2N}} \ \ \ \ \ \ (N>>1) \end{eqnarray}

La dispersión, para un caso típico de 100 mediciones es del orden de un 10% y con 1000 mediciones bajamos casi un orden de magnitud. Como regla del pistolero, con 100 mediciones podemos obtener dos cifras significativas de $\sigma$ y con 1000 mediciones tenemos una cifra más. En la práctica solemos reportar sólo una cifra significativa, porque lo único que nos interesa de esta dispersión es cómo truncar las estimaciones de $\mu$.

Puede sonar a poco, pero con estos resultado se vive, se come, se educa.

Hay que ver la distribución una vez en la vida:

In [4]:
def ch2i(sigma,N):
    '''Distribución chi^2 inversa
    '''
    if sigma==0: return 0
    
    lnout = -N*log(sigma) - N/2/sigma**2
    return exp(lnout)

figure(figsize=(10,6)) 
title('Probabilidad de $\sigma$, dados $N$ datos', fontsize=20)
xlabel('$\sigma/\sqrt{Var(x)}$', fontsize=20)
s = linspace(0.6,2,500)
out10 = array([ch2i(item,N=10) for item in s])
out10 = out10 /sum(out10)
plot(s,out10,label='$N=10$',lw=2,color='b')

out100 = array([ch2i(item,N=100) for item in s])
out100 = out100 /sum(out100)
plot(s,out100,label='$N=100$',lw=2,color='r')

out500 = array([ch2i(item,N=500) for item in s])
out500 = out500 /sum(out500)
plot(s,out500,label='$N=500$',lw=2,color='g')

xlim(0.6,2)
grid()
legend(fontsize=20);

Una reflexión final

Fíjense que puedo haber medido miles de datos, pero para estimar los parámetros de la distribución, me alcanza con conocer tres números: $$\overline{x}\ ,\ Var(x)\ ,\ N \ .$$ Se dice que estos números son una estadística suficiente en el contexto de este modelo.

2. Distribución de Poisson

Esta es una distribución de variable discreta, que toma valores sobre números enteros positivos.

Es la que describe el número de eventos que se observarán, si estos ocurren con un ritmo promedio conocido $\lambda$ y la ocurrencia de un evento no afecta la ocurrencia de otros eventos.

Más precisamente, si

  • $k$ es el número de veces que ocurre un evento y toma valores enteros positivos.
  • la ocurrencia de un evento no afecta la ocurrencia del evento siguiente.
  • la tasa a la ocurren los eventos es constante.
  • no puede haber dos eventos exactamente simultáneos.
  • la probabilidad en un intervalo es proporcional al largo del intervalo.

Se dice que una variable que satisface estas características es un proceso de Poisson. Hay un lindo artículo en la Wikipedia sobre este asunto.

Poisson la propuso para describir la llegada de expedientes por día en un juzgado, pero anda igual de bien para describir el número de llamadas telefónicas en un dado intervalo de tiempo, el número de mails recibidos por día, cuántos decaimientos nucleares mido por intervalo de tiempo, o cuántos eventos hay que generan picos en cierto rango de altura (este último ejemplo es el resultado de cualquier medición hecha con un multicanal).

La densidad de probabilidad tiene la forma:

$$ Poisson(k,\lambda) = \frac{\lambda^k\ e^{-\lambda}}{k!} $$

El valor medio de la distribución es $\lambda$ y la varianza es también $\lambda$.

2.1 Dibujo

Es fácil ver y no tan fácil demostrar que a medida que $\lambda$ aumenta, la distribución se aproxima a una gaussiana (pero siempre sobre valores discretos).

In [5]:
from scipy.stats import poisson

figure(figsize=(10,7))
k=arange(50)
title('Distribución de Poisson',fontsize=20)
xlabel('$k$',fontsize=20)
ylabel('$Po(x,\lambda)$',fontsize=20)

plot(k,poisson.pmf(k,mu=5),'bo',label=r'$\lambda=5$',lw=2)
plot(k,poisson.pmf(k,mu=5),'b-',lw=2, alpha=.3)

plot(k,poisson.pmf(k,mu=10),'ro',label=r'$\lambda=10$',lw=2)
plot(k,poisson.pmf(k,mu=10),'r-',lw=2, alpha=.3)

plot(k,poisson.pmf(k,mu=30),'go',label=r'$\lambda=30$',lw=2)
plot(k,poisson.pmf(k,mu=30),'g-',lw=2, alpha=.3)

grid()
legend(fontsize=20);

2.2 Estimación de Parámetro

Dado un conjunto de datos, tenemos que estimar un parámetro ($\lambda$). Como vimos en la introducción, con priors uniformes resulta que:

$$ P(\lambda \ | D I) \propto \prod_i Poisson(k_i ,\lambda) = \frac{\lambda^{N\overline{k}}\ e^{-N\lambda}}{\prod_i k_i!} \propto Poisson(N\overline{k}, N\lambda) $$

Esto me interesa como función de $\lambda$. Hay que tener cuidado que esta distribución sólo podemos llamarla de Poisson si la evaluamos sobre los vaores discretos del primer argumento. La distribución en función de $\lambda$ se llama distribución Gamma.

El máximo de esta distribución se obtiene de derivar la expresión e igualarla a cero. Resulta:

\begin{eqnarray} \hat{\lambda} &=& \overline{k} \\ \sigma_{\lambda} &=& \sqrt {\frac{\overline{k}}{N}} \end{eqnarray}

Es este caso, la estadística suficiente es $$ \overline{k} \ ,\ N\ . $$

3. Distribución Exponencial

Es una distribución continua con un solo parámetro y forma exponencial.

Describe los tiempos entre eventos en un proceso de Poisson. Tiene la propiedad de perder la memoria: si en un dado período no ocurrió un evento, la probabilidad de allí en adelante tiene la misma distribución que había al comienzo.

Es interesante que la moda (el valor del máximo de la función densidad de probabilidad) siempre es cero y el valor medio depende el parámetro.

La densidad de probabilidad tiene la forma:

$$ Exp(x, \tau) = \frac{e^{-t/\tau}}{\tau} $$

Con esta forma ni hace falta hacer un dibujito.

3.1 Estimación de Parámetro

Para estimar el parámetro, hacemos la productoria y resulta:

$$ P(\tau \ | D I) \propto \prod_i Exp(t_i ,\tau) = \frac {e^{-N\overline{t}/\tau}} {\tau^N} $$

Esta función es otra vez una distribución Gamma y se puede calcular el valor máximo y la varianza:

\begin{eqnarray} \hat{\tau} &=& \overline{t} \\ \sigma_{\tau} &=& \frac{\overline{t}}{\sqrt{N}} \end{eqnarray}

La estadística suficiente es $$ \overline{t} \ ,\ N\ . $$

Lo que tiene que quedar en una mente distraída

[Normal](https://en.wikipedia.org/wiki/Normal_distribution) [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution) [Exponencial](https://en.wikipedia.org/wiki/Exponential_distribution)
Suma de varias contribuciones aleatorias Número de eventos independientes que ocurren con una tasa de probabilidad constante en el tiempo Diferencia de tiempos en procesos de Poisson.
Pérdida de memoria
prámetros $\mu\ \ \ \ $ valor medio
$\sigma\ \ \ \ \ $ dispersión		 $\lambda\ \ \ \ $ tasa de eventos $\tau\ $ tiempo característico
PDF $$\frac{1}{\sqrt{2\pi \sigma^2}} \ e^{-\frac{1}{2\sigma^2} (x-\mu)^2}$$ $$\frac{\lambda^k\ e^{-\lambda}}{k!}$$ $$\frac{e^{-t/\tau}}{\tau}$$
pomedio $\mu$ $\lambda$ $\tau$
moda (max) $\mu$ $[\lambda]-1\ , [\lambda]\ $ $0$
dipersión $\sigma$ $\sqrt{\lambda}$ $\tau$
estimadores $\hat{\mu} = \overline{x}\ \ \ \ \ \ \ \ \ \sigma_\mu = \frac{\sigma}{\sqrt{N}}$
$\hat{\sigma} = \sqrt{Var(x)}\ \ \ \sigma_\sigma = \frac{1}{\sqrt{2N}}$ 		$\hat{\lambda} = \overline{k}\ \ \ \ \ \ \ \ \ \sigma_\lambda = \frac{\overline{k}}{\sqrt{N}}$ $\hat{\tau} = \overline{t}\ \ \ \ \ \ \ \ \ \sigma_\tau = \frac{\overline{t}}{\sqrt{N}}$

4. Distribución Normal Multivariada

Necesitamos introducir un ejemplo de distribución de más de una variable.

En muchas aplicaciones, aparecen cosas que se parecen a campanas gaussianas, así que es el primer ejemplo que tenemos que entender con cierto detalle.

Lo más sencillo sería multiplicar las distribuciones gaussianas de dos variables independientes, y así vamos a comenzar.

4.1 Variables independientes

Tenemos las variables independientes $x$ e $y$ y queremos expresar la distribución conjunta. Es sencillo, porque son independientes:

\begin{eqnarray} P(x,y\ | I) &=& P(x\ | I)\ P(y\ | I) \\ &=& G(x,\mu_x,\sigma_x) \ G(y,\mu_y,\sigma_y) \\ &=& \frac{1}{\sqrt{2\pi \sigma_x^2}} \ \frac{1}{\sqrt{2\pi \sigma_y^2}}\ e^{-\frac{1}{2\sigma_x^2} (x-\mu_x)^2-\frac{1}{2\sigma_y^2} (y-\mu_y)^2} \end{eqnarray}

a esta última expresión, la podemos pensar en términos de una matriz diagonal $\mathbf{\Sigma}$ definida como:

$$ \mathbf{\Sigma} = \begin{pmatrix} \sigma_x^2 & 0 \\ 0 & \sigma_y^2 \end{pmatrix} \\ P(x,y\ | I) = Multivar\ (x,y,\mathbf{\mu},\mathbf{\Sigma}) = \frac{1}{\left(\sqrt{2\pi}\right)^2 \ \sqrt{\left|\mathbf{\Sigma}\right|}} \ e^{-\frac{1}{2}\ (x-\mu_x\ ;\ y-\mu_y ) \cdot \mathbf{\Sigma}^{-1} \cdot \begin{pmatrix} x-\mu_x \\ y-\mu_y\end{pmatrix}} $$

y así queda una forma muy compacta de una distribución de dos variables. Esta distribución se llama normal multivariada y se puede generalizar fácilmente a muchas variables.

El hecho de que sean independientes está reflejado en que la distribución conjunta se construye multiplicando la distribución de cada variable por separado. Visualmente, esto se ve en que la forma resultante está alineada según los ejes.

In [50]:
from scipy.stats import norm
from mpl_toolkits.mplot3d.axes3d import Axes3D

#### Parámetros
mux, muy = 0, 1
sx ,  sy = 1, .5

#### Funciones
def Gx(x):
    return norm.pdf(x,loc=mux,scale=sx)

def Gy(y):
    return norm.pdf(y,loc=muy,scale=sy)

def Joint(x,y):
    return Gx(x) * Gy(y)

#### Límites
xm, xM = mux-2*sx, mux+2*sx
ym, yM = muy-2*sy, muy+2*sy
xm = min(xm-mux,ym-muy)
xM = max(xM-mux,yM-muy)
ym,yM = muy+(xm-mux),muy+(xM-mux)
zM = Joint(mux,muy)

#### Grilla
X = linspace(xm, xM, 50)
Y = linspace(ym, yM, 50)
X, Y = meshgrid(X, Y)
Z = Joint(X,Y)

#### Dibujo
fig = plt.figure(figsize=(10,8))
ax  = fig.gca(projection='3d')

surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=1)
cset = ax.contourf(X, Y, Z, zdir='z', offset=-zM*1.5, cmap=cm.coolwarm)

#### Vista
ax.set_xlabel('$x$',fontsize=20)
ax.set_xlim(xm, xM)
ax.set_ylabel('$y$',fontsize=20)
ax.set_ylim(ym, yM)
ax.set_zlabel('$P\ (x,y)$',fontsize=20)
ax.set_zlim(-zM*1.5, zM);

ax.view_init(elev=25, azim=-58)           # elevation and angle
ax.dist=11

4.2 Variables correlacionadas

El ejemplo anterior es fácilmente generalizable, si liberamos la restricción que la matriz $\mathbf{\Sigma}$ sea diagonal. Lo vamos a relajar, pero no tanto: vamos a necesitar que la matriz sea simétrica y (semi) definida positiva.

Eso nos permitirá escribir en el exponente expresiones de segundo orden en $x$ e $y$.

Se puede demostrar que los elementos de nuestra matriz representan la covariancia entre las variables:

$$ \mathbf{\Sigma}_{ij} = Cov(x_i,x_j) = \overline{(x_i - \overline{x_i} ) \, (x_j - \overline{x_j})} $$

La covarianza es una medida de cuánto cambian a la vez dos variables. Si las variables son independientes, la covarianza es nula. La inversa no es estrictamente cierta, pero para los casos con los que nos vamos a encontrar en el barrio, vamos a usarlo como si así lo fuera.

Los elementos de la diagonal, son las correlaciones de una variable consigo misma y la definición coincide con la varianza.

$$ \mathbf{\Sigma} = \begin{pmatrix} \sigma_x^2 & \sigma_x \sigma_y \rho \\ \sigma_x \sigma_y \rho & \sigma_y^2 \end{pmatrix} $$

Aquí escribimos una versión normalizada de la covariancia, que es el coeficiente de correlación.

$$ \rho = \frac{Cov(x_i,x_j)}{\sigma_i \sigma_j} $$
  • Gracias a la normalización, $\rho$ varía entre $-1$ y $+1$.
  • En los extremos, las variables están tan correlacionadas que podemos saber el valor de una de ellas, conociendo la otra.
  • Si $\rho=0$ las variables son independientes (salvo casos muy raros).
  • Si la correlación es positiva, ambas variables aumentan o disminuyen juntas, si la correlación es negativa, cuando una variable aumenta, la otra disminuye.

En el caso general, la distribución normal multivariada no tiene nada de raro, salvo que el máximo no se encuentra orientado según los ejes. Se puede demostrar fácilmente que basta una rotación en el espacio $x,y$ para tener nuevas variables independientes.

Para ilustración, nada mejor que un dibujito. Reproducimos el caso anterior pero con una correlación.

In [66]:
from scipy.stats import multivariate_normal

#### Parámetros
mux, muy = 0, 1
sx ,  sy = 1, .5
rho = -.8

#### Matriz de covariancia
sigma = array([[sx**2 , sx*sy*rho],
               [sx*sy*rho , sy**2]])

#### Funciones
def Joint2(x,y):
    return multivariate_normal.pdf( array([x,y]),
                                     mean=array([mux,muy]),
                                     cov=sigma)

#### Límites
xm, xM = mux-2*sx, mux+2*sx
ym, yM = muy-2*sy, muy+2*sy
xm = min(xm-mux,ym-muy)
xM = max(xM-mux,yM-muy)
ym,yM = muy+(xm-mux),muy+(xM-mux)
zM = Joint2(mux,muy)

#### Grilla
X = linspace(xm, xM, 50)
Y = linspace(ym, yM, 50)
X, Y = meshgrid(X, Y)
Z = array([Joint2(x,y) for x,y in zip(X.ravel(),Y.ravel()) ]).reshape(shape(X))
#### Dibujo
fig = plt.figure(figsize=(10,8))
ax  = fig.gca(projection='3d')

surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=1)
cset = ax.contourf(X, Y, Z, zdir='z', offset=-zM*1.5, cmap=cm.coolwarm)

#### Vista
ax.set_xlabel('$x$',fontsize=20)
ax.set_xlim(xm, xM)
ax.set_ylabel('$y$',fontsize=20)
ax.set_ylim(ym, yM)
ax.set_zlabel('$P\ (x,y)$',fontsize=20)
ax.set_zlim(-zM*1.5, zM);

ax.view_init(elev=25, azim=-50)           # elevation and angle
ax.dist=11

4.3 Aproximación del Hesiano

Hay muchas distribuciones de probabilidad que tienen un pico simétrico y marcado, que solemos aproximar por una multinomial multivariada cerca del máximo. Eso nos permite estraer conclusiones sobre las distribuciones finales.

Lo que haremos es aproximar la distribución que nos interesa,

$$ P (\theta) \propto \exp \left(-\frac{\chi^2_{(\theta)}}{2}\right) $$

desarrollando a segundo orden el exponente alrededor del máximo:

$$ \chi^2_{(\theta)} \approx \chi^2_{(\hat{\theta})} + \frac{1}{2} (\theta-\hat{\theta})^T\ \mathbf{H} \,(\theta-\hat{\theta}) $$

donde:

  • Los términos de primer orden en $\theta$ se fueron porque estamos alrededor del máximo
  • La matriz $\mathbf{H}$ es la matriz de las derivadas segundas de $\chi^2$ y en el ambente se la llama el Hessiano

Comparando esta aproximación con la expresión de la distribución multinomial, podemos identificar que la matriz de correlación es:

$$ \mathbf{\Sigma} = 2\ \mathbf{H}^{-1} $$

y de esta manera construimos una aproximación a la distribución cerca del máximo:

$$ P (\theta) \approx Multivar\ (\theta,\hat{\theta},\mathbf{\Sigma}) $$

Resultado muy conveniente, por cierto y que usamos para hacer el ajuste de una recta.

_Y eso es todo, amigos . . ._