Campos aleatorios I: simulación de campos gaussianos

Campos aleatorios (random fields) I¶

Autor: Eduardo Martín Calleja

Los campos aleatorios, y especialmente los denominados campos aleatorios gaussianos (GRF), aparecen muy frecuentemente en el análisis estadístico en cosmología. Por ejemplo, el fondo cósmico de microondas (CMB) es un campo gaussiano en la superficie esférica, así como lo es aproximadamente la distribución de galaxias a gran escala, la distribución de la densidad de materia, etc. Por esta razón abordaré en este y sucesivos posts el estudio práctico de este tipo de campos. En particular, el objeto de esta primera publicación será la definición del concepto de campo aleatorio y campo aleatorio gaussiano, así como proporcionar un primer método para simular realizaciones generando ejemplos de este tipo de campos.

De forma imprecisa, un campo aleatorio $X$ es un conjunto de variables aleatorias indexadas espacialmente. Por ejemplo la altura de los puntos en una región del oceano, la presión atmosférica en los puntos de un área, el número de insectos en cada arbol de una masa forestal, o el valor de la fluctuación del CMB en cada dirección de la esfera celeste. Habitualmente los valores de la variable correspondiente en un punto del espacio no serán independientes del valor en otros puntos próximos, sino que existirá una relación entre las variables que vendrá plasmada en la función de covarianza que veremos más adelante. Es justamente esta covarianza o relación espacial entre variables la que dota de interés y complejidad a este tipo de construcciones matemáticas.

Como de costumbre, este post ha sido escrito íntegramente con el Notebook de IPython

Importaciones y referencias¶

In [1]:

%matplotlib inline

Declaraciones de funciones LATEX

\DeclareMathOperator{\Cov}{Cov}

\DeclareMathOperator{\Var}{Var}

$$\DeclareMathOperator{\Cov}{Cov}$$$$\DeclareMathOperator{\Var}{Var}$$

In [2]:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

%load_ext version_information
%version_information numpy, matplotlib

Out[2]:

Software	Version
Python	2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython	3.1.0
OS	Linux 3.13.0 54 generic x86_64 with debian jessie sid
numpy	1.9.2
matplotlib	1.4.3
Wed Jun 17 23:42:29 2015 CEST

In [3]:

np.set_printoptions(precision=4)  # Mostrar con solo cuatro decimales en notación científica
np.set_printoptions(suppress=True) # Mostrar como cero los valores muy pequeños en notación científica

Referencias:

En la primera parte, definiciones, etc. seguiré el libro de Eric Vanmarcke: "Random Fields, analysis and synthesis"

Esta es una página de enlaces muy completa

Y de gran utilidad me han resultado las "Lecture Notes" sobre métodos numéricos para la generación de campos aleatorios gaussianos, que se pueden descargar en la página web de Catherine E. Powell

Definiciones¶

A continuación se va a plantear un resumen de definiciones y resultados básicos que, aunque un poco extenso, nos servirá de referencia para este post y posts sucesivos.

Un campo aleatorio escalar real $\lbrace X(\mathbf t) \: : \: \mathbf t \in T \rbrace$ es una colección de variables aleatorias con valores reales, definidas sobre un espacio de probabilidad dado, indexadas por puntos de coordenadas $\mathbf t = (t_1, \dots, t_n)$ en un espacio n-dimensional de "parámetros", el cual puede ser por ejemplo un subconjunto de un espacio cartesiano $\mathbb R^n$. Supondremos $n>1$, ya que cuando $n=1$ se tiene un proceso estocástico. Una realización del campo aleatorio es el resultado colectivo de la obtención de un resultado (observación) del experimento aleatorio $X(\mathbf t)$ en cada punto $\mathbf t$ del espacio de parámetros. Designaremos a una realización particular con la notación $\lbrace x(\mathbf t) \rbrace$. Al resultado de la observación en un punto concreto $\mathbf t$ le llamaremos estado del campo aleatorio en esa ubicación.

A menudo eliminaremos las llaves en la notación, por lo que la expresión $X(\mathbf t)$ se podrá interpretar, en función del contexto, como una única variable aleatoria en el punto $\mathbf t$ del espacio de parámetros, o bien como la familia de todas las v.a. al recorrer $\mathbf t$ el espacio de parámetros.

Supongamos que seleccionamos $M$ localizaciones diferentes $(\mathbf t_1, \dots, \mathbf t_M)$ del campo aleatorio. Esto nos lleva a considerar el vector aleatorio $\mathbf X $ formado por las v.a. $X(\mathbf t_1), \dots, X(\mathbf t_M)$. Definiremos la función de distribución acumulativa conjunta (joint cdf) como la función que hace corresponder a cada $M$-tupla $\mathbf x = (x_1, \dots, x_M)$ de números reales el valor:

$$F_{\mathbf X}(\mathbf x) \equiv F(\mathbf x) \equiv F(x_1, \dots, x_M) = P(X(\mathbf t_1) \leq x_1, \dots, X(\mathbf t_M) \leq x_M )$$

A partir de ella, si el campo aleatorio es continuo, podemos definir la función de densidad de probabilidad conjunta (join pdf) del vector aleatorio $\mathbf X $ como la función definida como:

$$f_{\mathbf X}(\mathbf x) \equiv f(\mathbf x) \equiv f(x_1, \dots, x_M) = \frac{\partial^M F(x_1, \dots, x_M)}{\partial x_1, \dots, \partial x_M}$$

En el caso de un campo aleatorio discreto, en el cual cada variable aleatoria $X(\mathbf t)$ puede tomar solo valores discretos $x_1, x_2, \dots $, utilizaremos en su lugar el concepto de función de masa de probabilidad conjunta, definida como:

$$p_{\mathbf X}(\mathbf x) \equiv p(\mathbf x) \equiv p(x_1, \dots, x_M) = P(X(\mathbf t_1) = x_1, \dots, X(\mathbf t_M) = x_M )$$

Las variables aleatorias $X(\mathbf t_1), \dots, X(\mathbf t_M)$ se dicen independientes si y solo si:

$$ f(x_1, \dots, x_M) = f(x_1) \cdots f(x_M) $$

Un campo aleatorio se dice homogeneo o estacionario si todas las funciones de distribución conjuntas permanecen sin cambio cuando se aplica una traslación en el espacio de parámetros a las localizaciones $\mathbf t_1, \mathbf t_2, \dots$.

Un campo aleatorio se define como isótropo cuando las funciones de distribución conjuntas son invariantes frente a una rotación en el espacio de parámetros.

Momentos¶

Comencemos con una variable aleatoria $X$ con una función de densidad de probabilidad $f(x)$. Sea $g(X)$ una función de dicha v.a. Definiremos la esperanza de $g(X)$ como:

$$E[g(X)] = \int_{-\infty}^\infty g(x) f(x) dx $$

Que en el caso de una variable discreta se reduce a:

$$E[g(X)] = \sum_{i} g(x_i) p(x_i)$$

En el caso particular de que $g(X) = X^n$ obtendremos los momentos (absolutos) de la v.a. $X$:

Si $n=1$ tenemos el primer momento: esperanza o media de $X$:

$$\mu_X \equiv \mu = E[X] = \int_{-\infty}^\infty x f(x) dx $$

Si $n=2$ tenemos el segundo momento (mean square), cuya raiz cuadrada es la media cuadrática (root mean square, RMS).

$$E[X^2] = (X_{\text{RMS}})^2 = \int_{-\infty}^\infty x^2 f(x) dx $$

Y así sucesivamente. De forma análoga podemos definir los momentos centrales que se obtienen tomando $g(X) = (X-\mu)^n$. El primer momento central es cero, y el segundo momento central es la varianza de $X$:

$$\sigma_X^2 \equiv \sigma^2 = \Var[X] = E[(X-\mu)^2] = E[X^2] - \mu^2 $$

Su raiz cuadrada es la desviación estándar $\sigma$.

En el caso de un campo aleatorio, basta sustituir $X$ por $X(\mathbf t)$ , es decir, por la v.a. en una localización particular, en las definiciones anteriores. De esta manera los momentos anteriores serán en general una función del parámetro $\mathbf t$. Así escribiremos por ejemplo $\mu(\mathbf t)$, $\sigma(\mathbf t)$, etc. En general estos momentos dependerán de la posición, pero en el caso en que el campo aleatorio sea homogéneo serán constantes, y se podrá escribir:

$$E[X(\mathbf t)] = \mu(\mathbf t) = \mu$$$$\Var[X(\mathbf t)] = \sigma^2(\mathbf t) = \sigma^2$$

Covarianza¶

Comencemos definiendo la esperanza de una función $g(X_1, X_2)$ de dos variables aleatorias $X_2$ y $X_2$, cuya función de densidad conjunta es $f(x_1, x_2)$ como:

$$E[g(X_1, X_2)] = \int_{-\infty}^\infty \int_{-\infty}^\infty g(x_1, x_2) f(x_1, x_2) dx_1 dx_2$$

La covarianza de $X_1$ y $X_2$ se define como la esperanza del producto de las desviaciones de sus medias respectivas:

$$\Cov[X_1, X_2] \equiv \sigma_{12} = E[(X_1 - \mu_1)(X_2 - \mu_2)] = E[X_1 X_2] - \mu_1 \mu_2$$

Es evidente que, para una v.a. $X$ se cumple que:

$$\Var[X] = \sigma_{XX} = \Cov[X, X]$$

Si dividimos la covarianza por el producto de las desviaciones estándar, obtenemos una magnitud sin dimensiones, el coeficiente de correlación entre $X_1$ y $X_2$:

$$\rho_{12} = \frac{\Cov[X_1, X_2]}{\sigma_1 \sigma_2} = \frac{\sigma_{12}}{\sigma_1 \sigma_2}$$

Si $\Cov(X_1, X_2) =0 $ diremos que $X_1, \, X_2$ no están correlacionadas. Si dos variables son independientes, no están correlacionadas, aunque la inversa no es cierta en general. Lo es en el caso en que las dos variables sigan conjuntamente una distribución gausiana bidimensional.

Matriz de covarianzas¶

Supongamos un vector aleatorio $\mathbf X$ formado por $M$ variables aleatorias $X_1, \dots, X_M$. Definiremos la matriz de covarianza $\Sigma$ como:

$$ \Sigma_\mathbf X = \begin{pmatrix} \sigma_{11} & \sigma_{12} & \cdots & \sigma_{1M} \\ \sigma_{21} & \sigma_{22} & \cdots & \sigma_{2M} \\ \vdots & \vdots & & \vdots \\ \sigma_{M1} & \sigma_{M2} & \cdots & \sigma_{MM} \end{pmatrix} $$

En la que obviamente la diagonal principal de la matriz está formada por las varianzas de cada una de las variables.

$\Sigma$ es una matriz simétrica. Además, en el caso más general, es definida semipositiva, lo que significa que:

$$\forall \mathbf a \: = \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_M \end{pmatrix} \: \in \mathbb R^M, \: \mathbf a^T \, \Sigma \, \mathbf a \geq 0$$

O lo que es lo mismo:

$$\sum_{i=1}^M \sum_{j=1}^M a_i a_j \, \sigma_{ij} \geq 0$$

No obstante, exigiremos como propiedad adicional que $\Sigma$ sea definida positiva, lo cual quiere decir que, $$ \forall \mathbf a \neq 0, \: \mathbf a^T \, \Sigma \, \mathbf a > 0$$

Un caso en que la matriz no será definida positiva es cuando uno de los vectores es una combinación lineal de los otros. Normalmente esta situación no se presentará en el caso de los campos aleatorios que vamos a estudiar.

Las propiedades de la matriz de covarianzas (o de cualquier matriz simétrica definida positiva) son:

1. Sus valores propios $\lambda_1, \dots, \lambda_M$ son todos reales positivos: $\lambda_i > 0 \quad \forall i$ (si fuera semipositiva definida algunos valores propios podrian ser nulos, aunque nunca negativos)

2. Si representamos el vector aleatorio $\mathbf X$ como una matriz columna:

$$\mathbf X = \begin{pmatrix} X_1 \\ X_2 \\ \vdots \\ X_M \end{pmatrix} $$

es posible encontrar una transformación lineal de matriz $\mathbf C$ que transforma $\mathbf X$ en un vector aleatorio:

$$\mathbf Z = \mathbf C \mathbf X$$

cuya matriz de covarianza es diagonal y está formada por los valores propios de $\Sigma_\mathbf X$:

Esta matriz $\mathbf C$ satisface las siguientes propiedades:

$$\mathbf C \mathbf C^T = \mathbf I$$$$\Sigma_{\mathbf Z} = \mathbf C \Sigma_{\mathbf X} \mathbf C^T = \begin{pmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & \lambda_M \end{pmatrix} $$

Las filas de $\mathbf C$ son los vectores propios normalizados (norma unidad) de la matriz $\Sigma_{\mathbf X}$

Otra propiedad, que utilizaremos más adelante, es que si la matriz de covarianza es simétrica definida positiva, puede ser factorizada como el producto de una matriz triangular y su transpuesta, $\Sigma = L L^T$, siendo $L$ una matriz triangular, en la qque los términos por encima de la diagonal son nulos. A este tipo de matrices se les llama matrices de Cholesky

La función de covarianza de un campo aleatorio¶

Para obtener dos variables aleatorias en un campo aleatorio $X$, sobre las cuales calcular su covarianza, consideraremos el valor del campo en dos localizaciones diferentes: $X_1 = X(\mathbf t_1), \: X_2 = X(\mathbf t_2)$. La covarianza resultante será un valor real que dependerá de los dos puntos elegidos. Así definimos la función de covarianza de $X$ como:

$$C(\mathbf t, \mathbf t') \equiv C_X(\mathbf t, \mathbf t') = \Cov[X(\mathbf t), X(\mathbf t')] = E[X(\mathbf t)X(\mathbf t')] - \mu(\mathbf t)\mu(\mathbf t')$$

Si el campo es homogéneo, la función de covarianza solo dependerá de las posiciones relativas de los puntos $\mathbf t, \mathbf t'$, por lo que, en este caso, abusando de la notación definiremos una nueva función de una sola variable:

$$C(\boldsymbol \tau) = C(\mathbf t , \mathbf t') \: \forall \mathbf t , \mathbf t': \mathbf t - \mathbf t' = \boldsymbol \tau$$

Si el campo es además isotrópo, la función de covarianza dependerá solo de la distancia entre las dos localizaciones, y podemos definir:

$$C(\tau) = C(|\boldsymbol \tau|) = C(|\mathbf t - \mathbf t'|)$$

Distribución gaussiana¶

La función de densidad de probabilidad de una variable aleatoria gaussiana viene dada por:

$$f(x | \mu, \sigma) = \frac{1}{\sigma \sqrt{2 \pi}} \exp \left(-\frac{1}{2} \left(\frac{x-\mu}{\sigma}\right)^2 \right)$$

Cumpliendose que:

$$\mu_X = E[X] = \mu$$$$\sigma^2_X = E[(X-\mu)^2] = \sigma^2$$

Esta definición se generaliza para definir la distribución gaussiana multivariada.

Diremos que un vector aleatorio de dimensión $M$: $\mathbf X=(X_1, \dots, X_M)^T$ sigue una distribución multivariada gaussiana si para cada vector $\mathbf a = (a_1, \cdots, a_M)^T \in \mathbb R^M$, la variable unidimensional

$$\sum_{i=1}^d a_i X_i$$

es gaussiana.

En tal caso, existe un vector de medias $\boldsymbol \mu \in \mathbb R^M$, con $\mu_i = E[X_i]$ y una matriz de covarianzas $M \times M $, $\Sigma$

$$\Sigma_{i,j} \equiv \sigma_{ij} = E[(X_i - \mu_i)(X_j - \mu_j)]$$

Tales que la función de densidad de $\mathbf X$ es:

$$ f_\mathbf X(\mathbf x) = f_\mathbf X(x_1,\dots,x_M) = \frac{1}{(2 \pi)^{M/2} |\Sigma|^{1/2}} \: \exp \left( - \frac{1}{2} (\mathbf x-\boldsymbol \mu)^T \Sigma^{-1} (\mathbf x-\boldsymbol \mu) \right)$$

Donde $|\Sigma|$ es el determinante de $\Sigma$, $\Sigma^{-1}$ es la matriz inversa de $\Sigma$, $\mathbf x = (x_1,\dots, x_M)^T$, $\boldsymbol \mu = (\mu_1, \dots, \mu_M)^T$.

Escribiremos $\mathbf X \sim N(\boldsymbol \mu, \Sigma)$

Adviertase que la expresión anterior solo tiene sentido si la matriz de covarianzas es definida positiva, por lo que su determinante será no nulo.

Campos aleatorios gaussianos (GRF)¶

Un campo aleatorio $\lbrace X(\mathbf t) \: : \: \mathbf t \in T \rbrace$ es gaussiano si y solo si:

$\forall M > 0, \: \forall \, \mathbf t_1, \dots, \mathbf t_n \in T$, el vector aleatorio $\mathbf X = (X_1(\mathbf t_1), \dots, X_M(\mathbf t_M))^T$ sigue una distribución multivariada gaussiana.

De modo que, para cada uno de estos vectores aleatorios, podemos escribir $\mathbf X \sim N(\boldsymbol \mu, \Sigma)$, donde:

$$\boldsymbol \mu = (\mu_1, \dots, \mu_M)^T = (E[X(\mathbf t_1)], \dots, E[X(\mathbf t_M)]^T $$$$\Sigma = \left( \sigma_{ij} \right) \quad \text{con} \quad \sigma_{ij} = \Cov[X(\mathbf t_i), X(\mathbf t_j)] $$

Es decir, todas estas distribuciones multivariadas, y por lo tanto las distribuciones de probabilidad asociadas al campo aleatorio $X$, están totalmente determinadas por las dos funciones siguientes:

$$\mu(\mathbf t) = E[X(\mathbf t)]$$$$C(\mathbf t, \mathbf t') = \Cov[X(\mathbf t), X(\mathbf t')] = E[(X(\mathbf t)-\mu(\mathbf(t))(X(\mathbf t')-\mu(\mathbf t'))]$$

Esta es la propiedad fundamental y principal razón de la utilidad de los campos aleatorios gaussianos, el poderse definir a partir de la función de las esperanzas en cada punto del espacio de parámetros, y la covarianza en cada dos puntos de este mismo espacio.

Una propiedad, que no probaremos, de los campos aleatorios gaussianos es que, aunque la ausencia de correlación no implica en general la independencia de las variables, en el caso de los GRF si es cierto. Es decir, si $C(\mathbf t, \mathbf t') = 0$, se cumple que $X(\mathbf t)$ y $X(\mathbf t')$ son independientes y su función de densidad conjunta puede escribirse como: $f(x_1, x_2) = f(x_1)f(x_2)$.

En el caso de GRF estacionarios u homogeneos, como la función $\mu(\mathbf t)$ es constante, se suele suponer que $\mu =0$, ya que siempre se puede reconstruir el caso general sumando una constante al campo.

Tambien en este caso, se suele escribir la función de covarianza como función de una sola variable, lo que, fijando un punto arbitrario $\mathbf t_0$ del espacio de parámetros, y suponiendo como hemos comentado que $\mu=0$ nos permite escribir:

$$C(\mathbf t) = \Cov[X(\mathbf t_0), X(\mathbf t_0 + \mathbf t)] = E[(X(\mathbf t_0)X(\mathbf t_0 + \mathbf t)]$$

Ejemplos de funciones de covarianza¶

Veamos a continuación algunos ejemplos de funciones de covarianza que se usan a menudo para modelizar campos aleatorios gaussianos, y para hacer simulaciones de los mismos.

Frecuentemente se parte de funciones de covarianza de la forma:

$$C(\mathbf t, \mathbf t') = \sigma^2 V(\mathbf t, \mathbf t')$$

Con la propiedad de que $V(\mathbf t, \mathbf t) = 1$, es decir, que dan lugar a campos aleatorios gaussianos de varianza constante, o sea, la v.a. del campo en cada punto tiene la misma varianza. Esto se cumple en particular siempre que el campo es homogeneo.

En los ejemplos que siguen a continuación se emplean también unas longitudes de correlación $l_i, \: i = 1,\dots,n$, o medidas de distancia, que nos proporcionan una medida de la distancia a la que deben estar los puntos en cada dirección, para que las variables correspondientes estén correlacionadas. Todos estos ejemplos expresan la idea de que cuando los puntos están próximos, las variables correspondientes están más correlacionadas.

Así tenemos la función de covarianza exponencial separable:

$$C(\mathbf t, \mathbf t') = \sigma^2 \exp \left(- \sum_{i=1}^n \frac{|t_i - t_i'|}{l_i} \right) = \sigma^2 \prod_{i=1}^n \exp \left(-\frac{|t_i - t_i'|}{l_i} \right), \quad l_i > 0$$

La función de covarianza exponencial (que sería apropiada para un campo homogeneo e isotrópico, ya que solo depende de la distancia entre los dos puntos considerados):

$$C(\mathbf t, \mathbf t') = \sigma^2 \exp \left( - \frac{\parallel \mathbf t - \mathbf t' \parallel}{l} \right), \quad l>0$$

Otro ejemplo sería la función de covarianza gaussiana, que en su forma más general presenta longitudes de correlación diferentes en cada dimensión:

$$C(\mathbf t, \mathbf t') = \sigma^2 \exp \left(- \sum_{i=1}^n \frac{(x_i - x_i')^2}{l^2_i} \right)$$

No obstante, se suele utilizar una versión en la que las longitudes de correlación son iguales:

$$C(\mathbf t, \mathbf t') = \sigma^2 \exp \left( - \frac{\parallel \mathbf t - \mathbf t' \parallel^2}{l^2} \right), \quad l>0$$

Simulación de campos gaussianos¶

Lo primero será considerar campos aleatorios discretos. Para ello elegiremos $M$ parámetros $\mathbf t_1, \dots, \mathbf t_M$, y el correspondiente vector aleatorio formado por las variables aleatorias del campo en cada uno de estos puntos:

$$ \mathbf X = \begin{pmatrix} X(\mathbf t_1) \\ X(\mathbf t_2) \\ \vdots \\ X(\mathbf t_M) \end{pmatrix} = \begin{pmatrix} X_1 \\ X_2 \\ \vdots\\ X_M \end{pmatrix} $$

Junto con el correspondiente vector de medias:

$$ \boldsymbol \mu = \begin{pmatrix} \mu_1 \\ \mu_2 \\ \vdots \\ \mu_M \end{pmatrix} = \begin{pmatrix} E[X(\mathbf t_1)] \\ E[X(\mathbf t_2)] \\ \vdots \\ E[X(\mathbf t_M)] \end{pmatrix} $$

Y la matriz de covarianzas, la cual podemos escribir en forma matricial como:

$$\Sigma = E[(\mathbf X - \boldsymbol \mu)(\mathbf X - \boldsymbol \mu)^T]$$

En efecto:

$$E [(\mathbf X - \boldsymbol \mu)(\mathbf X - \boldsymbol \mu)^T] = E \left[ \begin{pmatrix} X_1 - \mu_1\\ X_2 - \mu_2\\ \vdots\\ X_M - \mu_M \end{pmatrix} \begin{pmatrix} X_1 - \mu_1, X_2 - \mu_2, \dots, X_M - \mu_M \end{pmatrix} \right] $$

Es decir, obtendremos:

$$ \begin{pmatrix} E[(X_1-\mu_1)(X_1-\mu_1)] & \cdots E[(X_1-\mu_1)(X_M-\mu_M)] \\ \vdots & \vdots \\ E[(X_M-\mu_M)(X_1-\mu_1)] & \cdots E[(X_M-\mu_M)(X_M-\mu_M)] \end{pmatrix} $$

Vamos ahora a ver como se puede generar una realización de un campo gaussiano discreto dado un vector de medias $\boldsymbol \mu$ y una matriz de covarianzas definida positiva $\Sigma$.

Lo primero será considerar $M$ variables aleatorias independientes $Z_1, Z_2, \dots, Z_M$ normales $N(0, 1)$. Es decir, el vector:

$$\boldsymbol Z = \begin{pmatrix} Z_1 \\ Z_2 \\ \vdots \\ Z_M \end{pmatrix}$$

Será un vector aleatorio gaussiano $\boldsymbol Z \sim N(\boldsymbol 0,I)$, siendo $I$ la matriz identidad. En efecto, las covarianzas fuera de la diagonal serán cero, por la independencia entre las variables.

A continuación, por ser $\Sigma$ simétrica definida positiva, admite una factorización como producto de una matriz y su traspuesta:

$$\Sigma = L L^T$$

Esta factorización puede hacerse de diversas formas, pero la más habitual es la descomposición de Cholesky, en la que la matriz $L$ es triangular inferior.

El vector aleatorio $\boldsymbol N = L \boldsymbol Z$ cumple: $\boldsymbol N \sim \boldsymbol N(\boldsymbol 0, \Sigma)$, por lo cual el vector $\boldsymbol X = \boldsymbol \mu + \boldsymbol N = \boldsymbol \mu + L \boldsymbol Z \sim N(\boldsymbol \mu, \Sigma) $. A partir de una realización del vector $\boldsymbol Z$ podremos construir una realización del vector $\boldsymbol X$.

En efecto:

$$E[\boldsymbol N] = E[L \boldsymbol Z] = L E[\boldsymbol Z] = \boldsymbol 0 $$$$E[\boldsymbol N \boldsymbol N^T] = E[L \boldsymbol Z \boldsymbol Z^T L^T] = L E[\boldsymbol Z \boldsymbol Z^T] L^T = L I L^T = \Sigma$$

Por lo que en efecto:

$$\boldsymbol N \sim N(\boldsymbol 0, \Sigma)$$

Un ejemplo con la factorización de Cholesky¶

A continuación se va a ver paso a paso un ejemplo completo de generación de una realización de un campo gaussiano sobre una región del plano, partiendo de una matriz de covarianza definida previamente. Esta matriz de covarianza se va a generar utilizando una función de covarianza de tipo exponencial, definida anteriormente, aunque para el caso podria servir cualquier matriz simétrica definida positiva. Sin embargo, la elección de este tipo de función de covarianza hará la simulación más "realista".

Para ello vamos a subdividir el cuadrado $[0,1] \times [0,1]$ mediente una red de $M \times M$ puntos. De esta forma el problema se reduce a calcular una realización del campo en los nodos de la red. Como función de covarianza elegiremos el tipo de covarianza exponencial definido anteriormente.

In [4]:

M = 40 # Número de puntos en cada dimensión en una subdivisión de [0,1]x[0,1]

#Utilizaremos una función de covarianza exponencial
# C(t1, t2) = sigma2 * exp(- norm(t1 - t2)/l)
#de parámetros:
sigma2 = 1    # Varianza (constante)
l = 1         # Longitud de correlación

# Los nodos t_n del dominio D = [0,1]x[0,1] se supondrán ordenados en una disposición secuencial:
#[(0,0), (0,1/M),..., (0,(M-1)/M), (1,0),...,(1,(M-1)/M),...,((M-1)/M, 0),...,((M-1)/M,(M-1)/M)]

#Secuencia de coordenadas de nodos en el mismo orden que aparecerán en
#las filas y columnas de la matriz de covarianza Sigma
nodos = [(i/M,j/M) for i in range(M) for j in range(M)]

In [5]:

len(nodos)

Out[5]:

1600

A continuación generaremos la matriz de covarianza, como un array de Numpy con dimensiones $(M*M, M*M)$

In [6]:

x,y = zip(*nodos)
X1, X2 = np.meshgrid(x,x)
Y1, Y2 = np.meshgrid(y,y)
# Utilizaremos una función de covarianza exponencial
Sigma = sigma2 * np.exp(-np.sqrt((X1 - X2)**2 +(Y1 - Y2)**2)/l)

In [7]:

Sigma.shape

Out[7]:

(1600, 1600)

El siguiente paso será la factorización de Sigma. Se utilizará la descomposición de Choleski, en que $\Sigma = L L^T$, siendo $L$ una matriz triangular inferior.

In [8]:

L = np.linalg.cholesky(Sigma)

Generamos ahora el vector normal $Z \sim N(\boldsymbol 0, I)$, el cual, al transformarlo linealmente mediante la matriz $L$, nos proporcionará una realización de un campo gaussiano $N(\boldsymbol 0, \Sigma)$. Otras realizaciones se podrían obtener generando un nuevo vector $Z$.

In [9]:

np.random.seed(1)    # Fijamos una semilla para hacer reproducible el resultado
Z = np.random.randn(len(nodos))
N = np.dot(L, Z)

In [10]:

len(N)

Out[10]:

1600

Finalmente vamos a representar gráficamente el array resultante con los valores del campo. Para ello generaremos un array bidimensional para situar en cada punto de la región cuadrada [0,1] x [0,1] el valor del campo escalar que corresponda

In [11]:

N2D = np.zeros((M,M))
for i in xrange(len(nodos)):
    x, y = divmod(i,M)
    N2D[x,y] = N[i]

Representación grafica de la realización simulada del campo mediante un gráfico de tipo "heatmap":

In [12]:

im = plt.imshow(N2D)
plt.colorbar(im);

Observese que la representación anterior no es cartesiana. Está indexada por el número de fila y número de columna del array. Unos gráficos más depurados se mostrarán en la sección siguiente.

Si se quisiera generar otra realización bastaría repetir el procedimiento anterior cambiando la semilla para generar un vector aleatorio $\boldsymbol Z$ diferente del anterior.

Influencia de los parámetros de la función de covarianza¶

A continuación, generaremos unos gráficos para mostrar la influencia de los parámetros $\sigma^2$ y $l$ de la función de covarianza exponencial en el resultado

In [13]:

def GRF(M, sigma2, l):
    '''Genera una realización aleatoria de un campo aleatorio gaussiano
    de vector de medias = 0, en la región [0,1]x[0,1] del plano
    utilizando una función de covarianza exponencial
    de parámetros sigma2, l, es decir:
    C(t1, t2) = sigma2 * exp(- norm(t1 - t2)/l)
    sigma2 es la varianza del campo en cada punto.
    l es la longitud de correlación.
    M es el número de subdivisiones de cada eje.
    La función devuelve un array de dimensiones (M, M) con el valor
    del campo escalar en cada punto
    '''
    # Los nodos t_n del dominio D = [0,1]x[0,1] se supondrán ordenados secuencialmente:
    #[(0,0), (0,1/M),...,(0,(M-1)/M),(1/M,0),...,(1/M,(M-1)/M),...,((M-1)/M, 0),...,((M-1)/M,(M-1)/M)]
    #generaremos una lista con las coordenadas de los nodos en el mismo orden que aparecerán en
    #las filas y columnas de la matriz de covarianza
    nodos = [(i/M,j/M) for i in range(M) for j in range(M)]
    #Generación de la matriz de covarianza Sigma
    x,y = zip(*nodos)
    X1, X2 = np.meshgrid(x,x)
    Y1, Y2 = np.meshgrid(y,y)
    Sigma = sigma2 * np.exp(-np.sqrt((X1 - X2)**2 +(Y1 - Y2)**2)/l)
    #Ahora hacemos la descomposición de Cholesky de Sigma
    L = np.linalg.cholesky(Sigma)
    np.random.seed(1)    #Fijamos una semilla para hacer reproducible el resultado
    Z = np.random.randn(len(nodos))    # Vector de M*M variables independientes N(0,1)
    N = np.dot(L, Z) #Generamos un vector aleatorio gaussiano N(0, Sigma) de longitud M*M
    N2D = np.zeros((M,M))
    for i in xrange(len(nodos)):
        a, b = divmod(i,M)    #indices en el array N2D
        N2D[a,b] = N[i]
    return N2D

In [14]:

M = 50
positions = np.linspace(0,M, 6)
labels = np.linspace(0,1,6)
fig, ax = plt.subplots(2,2, figsize=(15, 12))
im = ax[0,0].imshow(GRF(M, 1, 1).T, origin='lower')
ax[0,0].set_title(r'GRF $\sigma^2 = 1$, $l=1$', fontsize=16)
fig.colorbar(im, ax = ax[0,0])

im = ax[0,1].imshow(GRF(M, 1, 0.1).T, origin='lower')
ax[0,1].set_title(r'GRF $\sigma^2 = 1$, $l=0.1$', fontsize=16)
fig.colorbar(im, ax = ax[0,1])

im = ax[1,0].imshow(GRF(M, 9, 1).T, origin='lower')
ax[1,0].set_title(r'GRF $\sigma^2 = 9$, $l=1$', fontsize=16)
fig.colorbar(im, ax = ax[1,0])

im = ax[1,1].imshow(GRF(M, 9, 0.1).T, origin='lower')
ax[1,1].set_title(r'GRF $\sigma^2 = 9$, $l=0.1$', fontsize=16)
fig.colorbar(im, ax = ax[1,1])

for i in [0,1]:
    for j in [0,1]:
        ax[i,j].set_xticks(positions)
        ax[i,j].set_xticklabels(labels)
        ax[i,j].set_yticks(positions)
        ax[i,j].set_yticklabels(labels)

En las figuras anteriores se ve que la longitud de correlación $l$ en la función de covarianza tiene por efecto modificar la "influencia" de cada punto en los puntos vecinos. Un valor pequeño de $l$ resultará en un campo aleatorio con un aspecto más fragmentado.

En cambio, el parámetro $\sigma^2$ afecta a la varianza de la variable aleatoria en cada punto, considerada individualmente. La correlación entre puntos es la misma, como se aprecia en los gráficos, ya que los dos mapas correspondientes al mismo valor de $\sigma^2$ son prácticamente idénticos. Sin embargo se aprecia en la barra lateral que el rango de variación del campo escalar es mayor cuanto mayor es el valor de la varianza. Concretamente una varianza $\sigma^2 = 9$ resulta en un rango de variación de la variable aproximadamente $\sqrt 9 = 3$ veces mayor que el rango de variación para $\sigma^2 = 1$.

Para concluir, haremos una representación tridimensional de la simulación realizada. A la vista del gráfico que se muestra a continuación se comprenderá que los campos aleatorios gaussianos se empleen a menudo en aplicaciones de geoestadística y gráficos de ordenador para modelizar terrenos

In [15]:

M = 30
nodos = [(i/M, j/M) for i in range(M) for j in range(M)]
x = [i/M for i in range(M)]
X, Y = np.meshgrid(x,x)

In [16]:

fig = plt.figure(figsize=(7,7))
ax1 = fig.add_subplot(111, projection='3d')
N = GRF(M, 1, 1)
ax1.plot_surface(X, Y, N, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0.5)
ax1.set_zlim(-5, 5)
ax1.set_title(r'GRF $\sigma^2 = 1$, $l=1$', fontsize=16);

Termino en este punto el artículo, que una vez más ha resultado un tanto extenso. Espero que les resulte interesante. ¡Hasta la próxima!