Distribución propuesta para la distribución normal generalizada

Question

Yo soy la modelización de la dispersión de la planta el uso generalizado de la distribución normal (entrada de la wikipedia), que tiene la función de densidad de probabilidad:

$$ \frac{b}{2a\Gamma(1/b)} e^{-(\frac{d}{a})^b} $$

where $d$ is distance travelled, $$ is a scale parameter, and $b$ is the shape parameter. Mean distance travelled is given by the standard deviation of this distribution:

$$ \sqrt{\frac{a^2 \Gamma(3/b)}{\Gamma(1/b)}} $$

This is convenient because it allows for an exponential shape when $b=1$, a Gaussian shape when $b=2$, and for a leptokurtic distribution when $b<1$. This distribution crops up regularly in the plant dispersal literature, even though it is pretty rare in general, and hence difficult to find information about.

The most interesting parameters are $b$ and mean dispersal distance.

I am trying to estimate $$ and $b$ using MCMC, but I am struggling to come up with an efficient way to sample proposal values. So far, I have used Metropolis-Hastings, and drawn from uniform distributions $0 < a < 400 $ and $ 0 < b<3$, and I get posterior mean dispersal distances of about 200-400 metres, which does make biological sense. However, convergence is really slow, and I am not convinced it is exploring the full parameter space.

Its tricky to come up with a better proposal distribution for $un$ and $b$, because they depend on one another, without having much meaning on their own. The mean dispersal distance does have a clear biological meaning, but a given mean dispersal distance could be explained by infinitely many combinations of $$ and $b$. As such $$ and $b$ are correlated in the posterior.

So far I have used Metropolis Hastings, but I am open to any other algorithm that would work here.

Question: Can anyone suggest a more efficient way to draw proposal values for $$ and $b$?

Edit: información Adicional sobre el sistema: estoy estudiando una población de plantas a lo largo de un valle. El objetivo es determinar la distribución de las distancias recorridas entre por el polen entre las plantas de los donantes y las plantas que polinizan. Los datos que tengo son:

1. La ubicación y el ADN para cada posible donante de polen
2. Semillas recolectadas a partir de una muestra de 60 materna plantas (es decir, el polen receptores) que han sido cultivados y genotipo.
3. La ubicación y el ADN de cada materna de la planta.

No sé la identidad de las plantas donantes, pero esto puede ser inferido a partir de los datos genéticos mediante la determinación de que los donantes son los padres de cada plántula. Digamos que esta información está contenida en una matriz de probabilidades de G con una fila para cada una de las crías y una columna para cada uno de los candidatos de los donantes, que da la probabilidad de cada uno de los candidatos a ser el padre de cada uno de los descendientes basado en los datos genéticos sólo. G toma alrededor de 3 segundos para calcular, y debe ser recalculado en cada iteración, lo que ralentiza las cosas considerablemente.

Ya que generalmente se espera que cerca de candidatos a los donantes a ser más propensos a ser padres, paternidad inferencia es más precisa si usted conjuntamente inferir la paternidad y la dispersión. La matriz D tiene las mismas dimensiones como G, y contiene las probabilidades de paternidad basadas únicamente en función de las distancias entre la madre y el candidato y algunos vector de parámetros. La multiplicación de los elementos de D y G da la probabilidad conjunta de la paternidad dado genéticos y datos espaciales. El producto de la multiplicación de los valores da la probabilidad de dispersión de la modelo.

Como se describió anteriormente he estado usando el GND para el modelo de dispersión. De hecho, yo en realidad utiliza una mezcla de TIERRA y una distribución uniforme para permitir la posibilidad de que muy distante de los candidatos tener una mayor probabilidad de paternidad debido a la casualidad (la genética es desordenado, que iba a inflar la aparente cola de la MASA, si se ignora. Por lo que la probabilidad de la distancia de dispersión $d$ es:

$$ c \Pr(d|a,b) + \frac{(1-c)}{N}$$

where $\Pr(d|a,b)$ is the probability of dispersal distance from the GND, N is the number of candidates, and $c$ ($0< c <1$) determina la cantidad de la contribución de la MASA que hace a la dispersión.

Por consiguiente, existen dos consideraciones adicionales que aumentan la carga computacional:

1. Distancia de dispersión no se conoce, pero se deben deducir en cada iteración, y la creación de G para hacer esto es caro.
2. Hay un tercer parámetro, $c$, para integrar más.

Por estas razones es que a mí me parecía estar cada vez un poco demasiado complejo para realizar la cuadrícula de interpolación, pero estoy feliz de estar convencido de lo contrario.

Ejemplo

Aquí es un ejemplo simplificado de un código python que he utilizado. He simplificado de estimación de la paternidad de los datos genéticos, ya que esto implicaría una gran cantidad de código, y la reemplazó con una matriz de valores entre 0 y 1.

En primer lugar, definir funciones para calcular la MASA:

import numpy as np
from scipy.special import gamma

def generalised_normal_PDF(x, a, b, gamma_b=None):
    """
    Calculate the PDF of the generalised normal distribution.

    Parameters
    ----------
    x: vector
        Vector of deviates from the mean.
    a: float
        Scale parameter.
    b: float
        Shape parameter
    gamma_b: float, optional
        To speed up calculations, values for Euler's gamma for 1/b
        can be calculated ahead of time and included as a vector.
    """
    xv = np.copy(x)
    if gamma_b:
        return (b/(2 * a * gamma_b ))      * np.exp(-(xv/a)**b)
    else:
        return (b/(2 * a * gamma(1.0/b) )) * np.exp(-(xv/a)**b)

def dispersal_GND(x, a, b, c):
    """
    Calculate a probability that each candidate is a sire
    assuming assuming he is either drawn at random form the
    population, or from a generalised normal function of his
    distance from each mother. The relative contribution of the
    two distributions is controlled by mixture parameter c.

    Parameters
    ----------
    x: vector
        Vector of deviates from the mean.
    a: float
        Scale parameter.
    b: float
        Shape parameter
    c: float between 0 and 1.
        The proportion of probability mass assigned to the
        generalised normal function.
    """    
    prob_GND = generalised_normal_PDF(x, a, b)
    prob_GND = prob_GND / prob_GND.sum(axis=1)[:, np.newaxis]

    prob_drawn = (prob_GND * c) + ((1-c) / x.shape[1])
    prob_drawn = np.log(prob_drawn)

    return prob_drawn

Siguiente simular 2000 candidatos, y 800 descendencia. También simular una lista de las distancias entre las madres de los hijos y el candidato de los padres, y un maniquí G de la matriz.

n_candidates = 2000 # Number of candidates in the population
n_offspring  = 800 # Number of offspring sampled.
# Create (log) matrix G.
# These are just random values between 0 and 1 as an example, but must be inferred in reality.
g_matrix  = np.random.uniform(0,1, size=n_candidates*n_offspring)
g_matrix  = g_matrix.reshape([n_offspring, n_candidates])
g_matrix  = np.log(g_matrix)
# simulate distances to ecah candidate father
distances = np.random.uniform(0,1000, 2000)[np.newaxis]

Conjunto inicial de valores de parámetro:

# number of iterations to run
niter= 100
# set intitial values for a, b, and c.
a_current = np.random.uniform(0.001,500, 1)
b_current = np.random.uniform(0.01,  3, 1)
c_current = np.random.uniform(0.001,  1, 1)
# set initial likelihood to a very small number
lik_current = -10e12

Actualización a, b y c, y a su vez, y calcular la Metrópoli relación.

# number of iterations to run
niter= 100
# set intitial values for a, b, and c.
# When values are very small, this can cause the Gamma function to break, so the limit is set to >0.
a_current = np.random.uniform(0.001,500, 1)
b_current = np.random.uniform(0.01,  3, 1)
c_current = np.random.uniform(0.001,  1, 1)
# set initial likelihood to a very small number
lik_current = -10e12 
# empty array to store parameters
store_params = np.zeros([niter, 3])

for i in range(niter):
    a_proposed = np.random.uniform(0.001,500, 1)
    b_proposed = np.random.uniform(0.01,3, 1)
    c_proposed = np.random.uniform(0.001,1, 1)

    # Update likelihood with new value for a
    prob_dispersal = dispersal_GND(distances, a=a_proposed, b=b_current, c=c_current)
    lik_proposed = (g_matrix + prob_dispersal).sum() # lg likelihood of the proposed value
    # Metropolis acceptance ration for a
    accept = bool(np.random.binomial(1, np.min([1, np.exp(lik_proposed - lik_current)])))
    if accept:
        a_current = a_proposed
        lik_current = lik_proposed
    store_params[i,0] = a_current

    # Update likelihood with new value for b
    prob_dispersal = dispersal_GND(distances, a=a_current, b=b_proposed, c=c_current)
    lik_proposed = (g_matrix + prob_dispersal).sum() # log likelihood of the proposed value
    # Metropolis acceptance ratio for b
    accept = bool(np.random.binomial(1, np.min([1, np.exp(lik_proposed - lik_current)])))
    if accept:
        b_current = b_proposed
        lik_current = lik_proposed
    store_params[i,1] = b_current

    # Update likelihood with new value for c
    prob_dispersal = dispersal_GND(distances, a=a_current, b=b_current, c=c_proposed)
    lik_proposed = (g_matrix + prob_dispersal).sum() # lg likelihood of the proposed value
    # Metropolis acceptance ratio for c
    accept = bool(np.random.binomial(1, np.min([1, np.exp(lik_proposed - lik_current)])))
    if accept:
        c_current = c_proposed
        lik_current = lik_proposed
    store_params[i,2] = c_current

Answer 1

Usted no necesita usar la Cadena de Markov Monte Carlo (MCMC) método.

Si usted está usando el uniforme antes de distribuciones entonces usted está haciendo algo muy similar como la estimación de máxima verosimilitud, en un espacio restringido para los parámetros de $a$ e $b$.

$$P(a,b;d) = P(d;a,b) \frac{P(a,b)}{P(d)} = \mathcal{L}(a,b;d) \times const $$

donde $\frac{P(a,b)}{P(d)}$ es una constante (independiente de $a$ e $b$) y que se pueden encontrar por la escala de la probabilidad de la función que se integra a 1.

El registro de probabilidad de la función, para $n$ iid variables $d_i \sim GN(0,a,b)$ es:

$$\log \mathcal{L}(a,b;d) = -n \log(2a) - n \log\left(\frac{\Gamma(1/b)}{b} \right) - \frac{1}{a^b} \sum_{i=1}^n \left( d_i \right)^b $$

Para esta función no debería ser demasiado difícil de trazar y encontrar un máximo.

Answer 2

No acabo de entender cómo se está configurando el modelo: en particular, me parece que para una determinada semilla, la posible dispersión polínica distancias son un conjunto finito, y por lo tanto su "dispersión de la probabilidad" podría ser mejor llamar a una "tasa de dispersión" (como lo sería necesario normalizar sumando más putativo a los padres a ser una probabilidad). Por lo tanto los parámetros no tienen el significado (como en, valores plausibles) que usted espera.

He trabajado en un par de problemas similares en el pasado y por lo que voy a tratar de llenar los vacíos en la comprensión, como una manera de lo que sugiere un posible enfoque o mirada crítica. Le pido disculpas si me se pierda por completo el punto de tu pregunta original. El tratamiento que a continuación sigue, básicamente, Hadfield et al (2006), uno de los mejores artículos acerca de este tipo de modelo.

Deje $X_{l,k}$ denotar el genotipo en el locus $l$ para algunas de las $k$. Para la descendencia $i$ conocidos con la madre $m_i$ , y su padre putativo $f$, dejar que la probabilidad de que el observado crías genotipos ser $$G_{i,f} = \prod_l \mathrm{Pr}(X_{l,i}|X_{l,m_i}, X_{l,f}, \theta)$$ -- in the simplest case this is simply a product of Mendelian inheritance probabilities, but in more complicated cases may include some model of genotyping error or missing parental genotypes, so I include nuisance parameter(s) $\theta$.

Deje $\delta_i$ ser el polen de la distancia de dispersión de la descendencia $i$, y deje $d_{m_i,f}$ ser la distancia entre la madre $m_i$ , y su padre putativo $f$, y deje $D_{i,f} = q(d_{m_i,f}|a,b,c)$ ser la tasa de dispersión (e.g es una combinación ponderada de la generalizada normal y uniforme de documentos pdf como en tu pregunta). Para expresar la tasa de dispersión como una probabilidad, normalizar w.r.t para el estado finito en el espacio: el (finito) conjunto de posibles distancias de dispersión inducida por el (finito) número de supuestos padres en su área de estudio, por lo que $$\tilde{D}_{i,f} = \mathrm{Pr}(\delta_i = d_{m_i,f}|a,b,c) = \frac{D_{i,f}}{\sum_k D_{i,k}}$$

Deje $P_i$ ser la paternal asignación de semilla $i$, que es $P_i = f$ si la planta de $f$ es el padre de descendencia $i$. Suponiendo un uniforme antes sobre la paternidad de las asignaciones, $$\mathrm{Pr}(P_i = f|a,b,c,\theta,X) = \frac{G_{i,f}\tilde{D}_{i,f}}{\sum_k G_{i,k}\tilde{D}_{i,k}} = \frac{G_{i,f}D_{i,f}}{\sum_k G_{i,k}D_{i,k}}$$ En otras palabras, condicional en otros parámetros y genotipos, paterna de asignación es un discreto r.v. con finito de apoyo, que se normaliza por la integración a través de dicho apoyo (posibles padres).

Así que una manera razonable para escribir un simple muestrario de este problema es la Metrópoli-dentro-de-Gibbs:

1. Condicional en $\{a,b,c,\theta\}$, la actualización de la paternidad asignaciones $P_i$ para todos los $i$. Este es un discreto r.v. con finito de apoyo de modo que usted puede fácilmente sacar una exacta de la muestra
2. Condicional en $\{P_i,\theta\}$, actualización de $a,b,c$ con una Metropolis-Hastings actualización. Para formar el destino, sólo el $D$ valores en las ecuaciones anteriores deben ser actualizado, así que esto no es costoso
3. Condicional en $\{P_i,a,b,c\}$, actualización de $\theta$ con un MH de actualización. Para formar el destino, el $G$ valores necesitan ser actualizados, lo cual es costoso, pero el $D$ no.

Para disminuir el costo de extraer muestras de $\{a,b,c\}$, puede realizar los pasos 1 y 2 varias veces antes de 3. Para afinar la propuesta de distribución en los pasos 2-3, usted podría utilizar muestras de un preliminar de ejecución para la estimación de la covarianza de la articulación posterior distribución por $\{a,b,c,\theta\}$. A continuación, utilice esta covarianza estimación dentro de un multivariante de Gauss propuesta. Estoy seguro de que este no es el enfoque más eficiente, pero es fácil de implementar.

Ahora, este esquema puede ser cerca de lo que ya lo están haciendo (yo no puedo decir cómo está la modelización de la paternidad de tu pregunta). Pero más allá computacional preocupaciones, mi punto más importante es que los parámetros de $a,b,c$ no puede tener el significado de que usted piensa que hacer, con respecto a la distancia media de dispersión. Esto es debido a que, en el contexto de la paternidad modelo de $\mathrm{Pr}(P_i|\cdot)$ he descrito anteriormente, $a,b,c$ entrar en tanto el numerador y el denominador (la normalización constante): por lo tanto, el arreglo espacial de las plantas tienen un potencial fuerte efecto en el que los valores de $a,b,c$ tienen una alta probabilidad o probabilidad posterior. Esto es especialmente cierto cuando la distribución espacial de las plantas es desigual.

Por último, le sugiero que eche un vistazo a ese Hadfield papel vinculado anteriormente y el acompañamiento de paquete de R ("MasterBayes"), si usted no tiene ya. Al menos puede proporcionar ideas.