10
$\begingroup$

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
$\endgroup$
9
  • 2
    $\begingroup$ Are you looking for a prior on a and b, or for a proposal distribution in a Metropolis-Hastings algorithm? You seem to have used both terms interchangeably. $\endgroup$ Sep 26, 2018 at 14:32
  • $\begingroup$ You are right - sorry for not being clear. I am most interested in a proposal distribution for MH. I have changed the title where I mentioned priors accordingly. $\endgroup$
    – tellis
    Sep 27, 2018 at 8:26
  • $\begingroup$ Under a flat or Jeffreys prior on $a$, i.e., $\pi(a)\propto 1$ or $\pi(a)\propto 1/a$ I believe that a change of variable to $\alpha=a^{-b}$ produces a closed-form conditional $\pi(a|b,\text{data})$. $\endgroup$
    – Xi'an
    Sep 27, 2018 at 19:36
  • $\begingroup$ See en.wikipedia.org/wiki/Generalized_normal_distribution $\endgroup$ Sep 27, 2018 at 20:07
  • $\begingroup$ It remains quite unclear whether or not you are interested in setting a prior that helps or in running Metropolis-Hastings more efficiently. $\endgroup$
    – Xi'an
    Oct 1, 2018 at 11:06

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ This gird interpolation would work for two parameters and observed distances, and might be what I end up doing. I am in fact doing joint estimation of dispersal distance and paternity inference, which involves at least one other parameter to integrate over, and an extra likelihood term which is really slow (~3 seconds per iteration) that really slows down the chain. I think I would need about 10-fold more iterations than I currently used for the markov chain. $\endgroup$
    – tellis
    Sep 28, 2018 at 12:10
  • $\begingroup$ @tellis those terms 'dispersal distance' and 'paternity inference' I do not really understand. Maybe you could provide some more concrete information by adding a dataset or a piece of it. While doing that you might as well speak more about the 'one other parameter'. So, what data is it that you do have? $\endgroup$ Sep 28, 2018 at 12:39
  • 1
    $\begingroup$ I have added an example using simulated data. $\endgroup$
    – tellis
    Oct 2, 2018 at 7:34
0
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ My approach is indeed modelled on Hadfield's, with two major changes: (1) seeds from a mother may be full siblings, and therefore not independent. The problem is therefore one of jointly inferring dispersal, paternity, bat also sibship structure as well. (2) I am using a fractional paternity approach to consider all candidates simultaneously in proportion to their likelihood of paternity, rather than update paternity assignments sequentially, because there is a large space of possible fathers to explore. $\endgroup$
    – tellis
    Oct 4, 2018 at 4:32
  • $\begingroup$ I am using the package FAPS to do those things. $\endgroup$
    – tellis
    Oct 4, 2018 at 4:33
  • $\begingroup$ My question is essentially asking about an efficient proposal distribution for do your point 2. The rest of your answer describes something very close to what I have already done, including the formulation of the product of G and D (but thanks for this - I wasn't sure I had done it correctly, so it is useful to know a second pair of eyes agree!). $\endgroup$
    – tellis
    Oct 4, 2018 at 4:49
  • $\begingroup$ I don't have a canned solution w.r.t proposal distribution, sorry. But I have a few observations: (1) Steps 1-2 are very cheap, and can be iterated many times with little cost before moving to step 3. Even with a shoddy proposal in step 2, a lot of iterations should suffice to "make large moves" in the state space of $a,b,c$. $\endgroup$
    – Nate Pope
    Oct 4, 2018 at 18:45
  • $\begingroup$ (2) The conditional distribution in step 2 is 3-dimensional. As in: easy to visualize. What does the unnormalized target of $a,b,c$ look like at a MAP estimate of the paternity assignments for a fixed $G$? Visualizing the unnormalized target across different paternities should give you a sense of if it's multimodal, flat in areas, etc. $\endgroup$
    – Nate Pope
    Oct 4, 2018 at 18:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.