I am modelling plant dispersal using a generalised normal distribution (wikipedia entry), which has the probability density function:
$$ \frac{b}{2a\Gamma(1/b)} e^{-(\frac{d}{a})^b} $$
where $d$ is distance travelled, $a$ 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 $a$ 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 $a$ 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 $a$ and $b$. As such $a$ 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 $a$ and $b$?
Edit: Additional information about the system: I'm studying a population of plants along a valley. The aim is to determine the distribution of distances travelled between by pollen between donor plants and the plants they pollinate. The data I have are:
- The location and DNA for every possible pollen donor
- Seeds collected from a sample of 60 maternal plants (i.e. pollen receivers) that have been grown and genotyped.
- The location and DNA for each maternal plant.
I do not know the identity of the donor plants, but this can be inferred from the genetic data by determining which donors are the fathers of each seedling. Let's say this information is contained in a matrix of probabilities G with a row for each offspring and a column for each candidate donor, that gives the probability of each candidate being the father of each offspring based on genetic data only. G takes about 3 seconds to compute, and needs to be recalculated at every iteration, which slows things down considerably.
Since we generally expect closer candidate donors to be more likely to be fathers, paternity inference is more accurate if you jointly infer paternity and dispersal. Matrix D has the same dimensions as G, and contains probabilities of paternity based only a function of distances between mother and candidate and some vector of parameters. Multiplying elements in D and G gives the joint probability of paternity given genetic and spatial data. The product of multiplied values gives the likelihood of dispersal model.
As described above I have been using the GND to model dispersal. In fact I actually used a mixture of a GND and a uniform distribution to allow for the possibility of very distant candidates having a higher likelihood of paternity due to chance alone (genetics is messy) which would inflate the apparent tail of the GND if ignored. So the probability of dispersal distance $d$ is:
$$ 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$) determines how much contribution the GND makes to dispersal.
There are therefore two additional considerations that increase computational burden:
- Dispersal distance is not known but must be inferred at each iteration, and creating G to do this is expensive.
- There is a third parameter, $c$, to integrate over.
For these reasons it seemed to me to be ever so slightly too complex to perform grid interpolation, but I am happy to be convinced otherwise.
Example
Here is a simplified example of the python code I have used. I have simplified estimation of paternity from genetic data, since this would involve a lot of extra code, and replaced it with a matrix of values between 0 and 1.
First, define functions to calculate the 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
Next simulate 2000 candidates, and 800 offspring. Also simulate a list of distances between the mothers of the offspring and the candidate fathers, and a dummy G matrix.
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]
Set initial parameter values:
# 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
Update a, b, and c in turn, and compute the Metropolis ratio.
# 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