Proposal distribution for a generalised normal distribution

Question

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:

1. The location and DNA for every possible pollen donor
2. Seeds collected from a sample of 60 maternal plants (i.e. pollen receivers) that have been grown and genotyped.
3. 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:

1. Dispersal distance is not known but must be inferred at each iteration, and creating G to do this is expensive.
2. 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

Answer 1

You do not need to use the Markov Chain Monte Carlo (MCMC) method.

If you are using uniform prior distributions then you are doing something very similar as maximum likelihood estimation on a restricted space for the parameters $a$ and $b$.

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

where $\frac{P(a,b)}{P(d)}$ is a constant (independent from $a$ and $b$) and it can be found by scaling the likelihood function such that it integrates to 1.

The log likelihood function, for $n$ iid variables $d_i \sim GN(0,a,b)$ is:

$$\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 $$

For this function it should not be too difficult to plot it and find a maximum.

Answer 2

I don't quite understand how you are setting up the model: in particular, it seems to me that for a given seed, the possible pollen dispersal distances are a finite set, and thus your "dispersal probability" might be better termed a "dispersal rate" (as it would need to be normalized by summing over putative fathers to be a probability). Thus the parameters may not quite have the meaning (as in, plausible values) that you expect.

I have worked on a couple of similar problems in the past and so I'll try to fill in the gaps in my understanding, as a way of suggesting a possible approach/critical look. Apologies if I completely miss the point your original question. The treatment below basically follows Hadfield et al (2006), one of the better papers about of this kind of model.

Let $X_{l,k}$ denote the genotype at locus $l$ for some individual $k$. For offspring $i$ with known mother $m_i$ and putative father $f$, let the probability of the observed offspring genotypes be $$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$.

Let $\delta_i$ be the pollen dispersal distance for offspring $i$, and let $d_{m_i,f}$ be the distance between known mother $m_i$ and putative father $f$, and let $D_{i,f} = q(d_{m_i,f}|a,b,c)$ be the dispersal rate (e.g a weighted combination of generalized normal and uniform pdfs as in your question). To express the dispersal rate as a probability, normalize w.r.t to the finite state space: the (finite) set of possible dispersal distances induced by the (finite) number of putative fathers in your study area, so that $$\tilde{D}_{i,f} = \mathrm{Pr}(\delta_i = d_{m_i,f}|a,b,c) = \frac{D_{i,f}}{\sum_k D_{i,k}}$$

Let $P_i$ be the paternal assignment of seed $i$, that is $P_i = f$ if plant $f$ is the father of offspring $i$. Assuming a uniform prior on paternity assignments, $$\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}}$$ In other words, conditional on other parameters and genotypes, paternal assignment is a discrete r.v. with finite support, that is normalized by integrating across said support (possible fathers).

So a reasonable way to write a simple sampler for this problem is Metropolis-within-Gibbs:

1. Conditional on $\{a,b,c,\theta\}$, update paternity assignments $P_i$ for all $i$. This is a discrete r.v. with finite support so you can easily draw an exact sample
2. Conditional on $\{P_i,\theta\}$, update $a,b,c$ with a Metropolis-Hastings update. To form the target, only the $D$ values in the equations above need to be updated, so this isn't costly
3. Conditional on $\{P_i,a,b,c\}$, update $\theta$ with a MH update. To form the target, the $G$ values need to be updated, which is costly, but the $D$ do not.

To decrease the cost of drawing samples of $\{a,b,c\}$, you could perform steps 1-2 multiple times before 3. To tune the proposal distributions in steps 2-3, you could use samples from a preliminary run to estimate the covariance of the joint posterior distribution for $\{a,b,c,\theta\}$. Then use this covariance estimate within a multivariate Gaussian proposal. I'm sure this isn't the most efficient approach, but it is easy to implement.

Now, this scheme may be close to what you are already doing (I can't tell how you are modelling paternity from your question). But beyond computational concerns, my larger point is that the parameters $a,b,c$ may not have the meaning you think they do, with regards to mean dispersal distance. This is because, in the context of the paternity model $\mathrm{Pr}(P_i|\cdot)$ I described above, $a,b,c$ enter into both numerator and denominator (normalizing constant): thus, the spatial arrangement of plants will have a potentially strong effect on which values of $a,b,c$ have a high likelihood or posterior probability. This is especially true when the spatial distribution of the plants is uneven.

Finally, I suggest you take a look at that Hadfield paper linked to above and the accompanying R package ("MasterBayes"), if you haven't already. At the least it may provide ideas.