Distribución de propuestas para una distribución normal generalizada.

Estoy modelando la dispersión de plantas usando una distribución normal generalizada ( entrada de wikipedia ), que tiene la función de densidad de probabilidad:

donde $d$ es la distancia recorrida, $a$ es un parámetro de escala y $b$ es el parámetro de forma. La distancia media recorrida viene dada por la desviación estándar de esta distribución:

Esto es conveniente porque permite una forma exponencial cuando $b=1$ , una forma gaussiana cuando $b=2$ y una distribución leptokurtica cuando $b<1$ . Esta distribución aparece regularmente en la literatura sobre dispersión de plantas, aunque es bastante rara en general y, por lo tanto, difícil de encontrar información.

Los parámetros más interesantes son $b$ distancia de dispersión media.

Estoy tratando de estimar $a$ y $b$ utilizando MCMC, pero estoy luchando para llegar a una forma eficiente de los valores de propuesta de la muestra. Hasta ahora, he usado Metropolis-Hastings, y extraído de distribuciones uniformes $0 < a < 400$ y $0 < b<3$ , y obtengo distancias de dispersión medias posteriores de aproximadamente 200-400 metros, lo que tiene sentido biológico. Sin embargo, la convergencia es realmente lenta, y no estoy convencido de que esté explorando el espacio de parámetros completo.

Su difícil de llegar a una mejor distribución propuesta de $a$ y $b$ , ya que dependen el uno del otro, sin tener mucho significado por sí mismos. La distancia de dispersión media tiene un significado biológico claro, pero una distancia de dispersión media dada podría explicarse por infinitas combinaciones de $a$ y $b$ . Como tal $a$ y $b$ están correlacionados en la parte posterior.

Hasta ahora he usado Metropolis Hastings, pero estoy abierto a cualquier otro algoritmo que funcione aquí.

Pregunta: ¿Puede alguien sugerir una manera más eficiente para extraer los valores de propuesta para $a$ y $b$ ?

Editar: 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 por el polen entre las plantas donantes y las plantas que polinizan. Los datos que tengo son:

1. La ubicación y el ADN de cada posible donante de polen
2. Semillas recolectadas de una muestra de 60 plantas maternas (es decir, receptores de polen) que han sido cultivadas y genotipadas.
3. La ubicación y el ADN de cada planta materna.

No conozco la identidad de las plantas donantes, pero esto se puede inferir de los datos genéticos al determinar qué donantes son los padres de cada plántula. Digamos que esta información está contenida en una matriz de probabilidades G con una fila para cada descendencia y una columna para cada donante candidato, que da la probabilidad de que cada candidato sea el padre de cada descendencia basándose solo en datos genéticos. G tarda unos 3 segundos en calcular, y necesita ser recalculado en cada iteración, lo que ralentiza considerablemente las cosas.

Como generalmente esperamos que los donantes candidatos más cercanos sean más propensos a ser padres, la inferencia de paternidad es más precisa si deduce conjuntamente la paternidad y la dispersión. La matriz D tiene las mismas dimensiones que G y contiene probabilidades de paternidad basadas solo en una función de distancias entre la madre y el candidato y algún vector de parámetros. Multiplicar elementos en D y G da la probabilidad conjunta de paternidad dada la información genética y espacial. El producto de valores multiplicados da la probabilidad de un modelo de dispersión.

Como se describió anteriormente, he estado usando GND para modelar la dispersión. De hecho, en realidad utilicé una mezcla de GND y una distribución uniforme para permitir la posibilidad de que candidatos muy distantes tengan una mayor probabilidad de paternidad debido solo al azar (la genética es desordenada), lo que inflaría la cola aparente de GND si se ignora. Entonces la probabilidad de distancia de dispersión $d$ es:

donde $\Pr(d|a,b)$ es la probabilidad de distancia de dispersión del GND, N es el número de candidatos $c$ ( $0< c <1$ ) determina cuánta contribución hace el GND a la dispersión.

Por lo tanto, hay dos consideraciones adicionales que aumentan la carga computacional:

1. La distancia de dispersión no se conoce, pero debe inferirse en cada iteración y crear G para hacer esto es costoso.
2. Hay un tercer parámetro, $c$ , para integrar.

Por estas razones, me pareció un poco demasiado complejo para realizar la interpolación de cuadrícula, pero estoy feliz de estar convencido de lo contrario.

Ejemplo

Aquí hay un ejemplo simplificado del código de Python que he usado. He simplificado la estimación de la paternidad a partir de datos genéticos, ya que esto implicaría una gran cantidad de código adicional, y lo reemplacé con una matriz de valores entre 0 y 1.

Primero, defina funciones para calcular el GND:

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

A continuación, simule 2000 candidatos y 800 descendientes. Simule también una lista de distancias entre las madres de la descendencia y los padres candidatos, y una matriz G simulada .

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]

Establecer valores de parámetros iniciales:

# 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

Actualice a, byc a su vez, y calcule la relación de Metrópolis.

# 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

No necesita utilizar el método Markov Chain Monte Carlo (MCMC).

$a$$b$

$\frac{P(a,b)}{P(d)}$$a$$b$

$n$$d_i \sim GN(0,a,b)$

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

No entiendo muy bien cómo está configurando el modelo: en particular, me parece que para una semilla determinada, las posibles distancias de dispersión del polen son un conjunto finito y, por lo tanto, su "probabilidad de dispersión" podría denominarse mejor " tasa de dispersión "(ya que tendría que normalizarse sumando padres supuestos para ser una probabilidad). Por lo tanto, los parámetros pueden no tener el significado (como en valores plausibles) que usted espera.

He trabajado en un par de problemas similares en el pasado, por lo que trataré de llenar los vacíos en mi entendimiento, como una forma de sugerir un posible enfoque / mirada crítica. Disculpas si pierdo completamente el punto de tu pregunta original. El tratamiento a continuación básicamente sigue a Hadfield et al (2006) , uno de los mejores artículos sobre este tipo de modelo.

$X_{l,k}$$l$$k$$i$$m_i$$f$

$\delta_i$$i$$d_{m_i,f}$$m_i$$f$$D_{i,f} = q(d_{m_i,f}|a,b,c)$

$P_i$$i$$P_i = f$$f$$i$

Entonces, una forma razonable de escribir una muestra simple para este problema es Metropolis-inside-Gibbs:

1. $\{a,b,c,\theta\}$$P_i$$i$
2. $\{P_i,\theta\}$$a,b,c$$D$
3. $\{P_i,a,b,c\}$$\theta$$G$$D$

$\{a,b,c\}$$\{a,b,c,\theta\}$

$a,b,c$$\mathrm{Pr}(P_i|\cdot)$$a,b,c$$a,b,c$tienen una alta probabilidad o probabilidad posterior. Esto es especialmente cierto cuando la distribución espacial de las plantas es desigual.

Finalmente, le sugiero que eche un vistazo al documento de Hadfield vinculado anteriormente y al paquete R que lo acompaña ("MasterBayes"), si aún no lo ha hecho. Al menos puede aportar ideas.