Angebotsverteilung für eine verallgemeinerte Normalverteilung

Ich modelliere die Ausbreitung von Pflanzen mithilfe einer verallgemeinerten Normalverteilung ( Wikipedia-Eintrag ), die die Wahrscheinlichkeitsdichtefunktion hat:

Dabei ist die zurückgelegte Strecke, ein Skalierungsparameter und der Formparameter. Die mittlere zurückgelegte Strecke ergibt sich aus der Standardabweichung dieser Verteilung:$d$$a$$b$

Dies ist praktisch, da es eine Exponentialform bei , eine Gaußsche Form bei und eine leptokurtische Verteilung bei . Diese Verteilung taucht regelmäßig in der Literatur zur Pflanzenverbreitung auf, obwohl sie im Allgemeinen ziemlich selten ist und daher schwer zu finden ist.$b=1$$b=2$$b<1$

Die interessantesten Parameter sind und die mittlere Ausbreitungsentfernung.$b$

Ich versuche, und mithilfe von MCMC zu schätzen , habe jedoch Schwierigkeiten, eine effiziente Methode zum Abtasten von Angebotswerten zu finden. Bisher habe ich Metropolis-Hastings verwendet und aus gleichmäßigen Verteilungen und , und ich erhalte hintere mittlere Ausbreitungsentfernungen von etwa 200 bis 400 Metern, was biologisch sinnvoll ist. Die Konvergenz ist jedoch sehr langsam, und ich bin nicht davon überzeugt, dass sie den gesamten Parameterraum untersucht.$a$$b$$0 < a < 400$$0 < b<3$

Es ist schwierig, eine bessere Angebotsverteilung für und , da sie voneinander abhängig sind, ohne selbst viel Bedeutung zu haben. Die mittlere Ausbreitungsentfernung hat eine klare biologische Bedeutung, aber eine gegebene mittlere Ausbreitungsentfernung könnte durch unendlich viele Kombinationen von und . Als solche sind und im posterioren Bereich korreliert.$a$$b$b a b$a$$b$$a$$b$

Bisher habe ich Metropolis Hastings verwendet, aber ich bin offen für jeden anderen Algorithmus, der hier funktionieren würde.

Frage: Kann jemand einen effizienteren Weg vorschlagen, um Vorschlagswerte für und zu zeichnen ?$a$$b$

Bearbeiten: Zusätzliche Informationen zum System: Ich untersuche eine Pflanzenpopulation entlang eines Tals. Ziel ist es, die Verteilung der Entfernungen zwischen Pollen zwischen Spenderpflanzen und den Pflanzen, die sie bestäuben, zu bestimmen. Die Daten, die ich habe, sind:

1. Der Ort und die DNA für jeden möglichen Pollenspender
2. Samen, die aus einer Probe von 60 Mutterpflanzen (dh Pollenempfängern) entnommen wurden, die gezüchtet und genotypisiert wurden.
3. Der Ort und die DNA für jede mütterliche Pflanze.

Ich kenne die Identität der Spenderpflanzen nicht, aber dies kann aus den genetischen Daten abgeleitet werden, indem bestimmt wird, welche Spender die Väter jedes Sämlings sind. Angenommen, diese Informationen sind in einer Wahrscheinlichkeitsmatrix G mit einer Zeile für jeden Nachwuchs und einer Spalte für jeden Spenderkandidaten enthalten, die die Wahrscheinlichkeit angibt, dass jeder Kandidat der Vater jedes Nachwuchses ist, nur basierend auf genetischen Daten. Die Berechnung von G dauert ungefähr 3 Sekunden und muss bei jeder Iteration neu berechnet werden, was die Dinge erheblich verlangsamt.

Da wir im Allgemeinen davon ausgehen, dass engere Spenderkandidaten eher Väter sind, ist die Schlussfolgerung der Vaterschaft genauer, wenn Sie gemeinsam auf Vaterschaft und Streuung schließen. Matrix D hat die gleichen Dimensionen wie G und enthält Vaterschaftswahrscheinlichkeiten, die nur auf einer Funktion der Abstände zwischen Mutter und Kandidat und einigen Parametervektoren beruhen. Die Multiplikation von Elementen in D und G ergibt die gemeinsame Wahrscheinlichkeit einer Vaterschaft bei genetischen und räumlichen Daten. Das Produkt multiplizierter Werte gibt die Wahrscheinlichkeit eines Ausbreitungsmodells an.

Wie oben beschrieben, habe ich die GND verwendet, um die Ausbreitung zu modellieren. Tatsächlich habe ich eine Mischung aus einer GND und einer gleichmäßigen Verteilung verwendet, um die Möglichkeit zu ermöglichen, dass sehr weit entfernte Kandidaten eine höhere Wahrscheinlichkeit der Vaterschaft haben, allein aufgrund des Zufalls (die Genetik ist chaotisch), was den scheinbaren Schwanz der GND aufblasen würde, wenn sie ignoriert würden. Die Wahrscheinlichkeit der Ausbreitungsentfernung ist also:$d$

Dabei ist die Wahrscheinlichkeit der Ausbreitungsentfernung von der GND, N die Anzahl der Kandidaten und ( ) bestimmt, wie viel Beitrag die GND zur Ausbreitung leistet.$\Pr(d|a,b)$$c$$0< c <1$

Es gibt daher zwei zusätzliche Überlegungen, die den Rechenaufwand erhöhen:

1. Die Ausbreitungsentfernung ist nicht bekannt, muss jedoch bei jeder Iteration abgeleitet werden, und das Erstellen von G , um dies zu tun, ist teuer.
2. Es gibt einen dritten Parameter, , über den integriert werden soll.$c$

Aus diesen Gründen schien es mir etwas zu komplex, um eine Gitterinterpolation durchzuführen, aber ich bin froh, davon überzeugt zu sein.

Beispiel

Hier ist ein vereinfachtes Beispiel für den von mir verwendeten Python-Code. Ich habe die Schätzung der Vaterschaft anhand genetischer Daten vereinfacht, da dies viel zusätzlichen Code beinhalten würde, und ihn durch eine Wertematrix zwischen 0 und 1 ersetzt.

Definieren Sie zunächst Funktionen zur Berechnung des 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

Als nächstes simulieren Sie 2000 Kandidaten und 800 Nachkommen. Simulieren Sie auch eine Liste der Abstände zwischen den Müttern der Nachkommen und den Kandidatenvätern sowie eine 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]

Anfängliche Parameterwerte einstellen:

# 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

Aktualisieren Sie nacheinander a, b und c und berechnen Sie das Metropolis-Verhältnis.

# 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

Sie müssen die Markov Chain Monte Carlo (MCMC) -Methode nicht verwenden.

Wenn Sie einheitliche vorherige Verteilungen verwenden, führen Sie eine sehr ähnliche Schätzung der maximalen Wahrscheinlichkeit auf einem begrenzten Raum für die Parameter und .$a$$b$

Dabei ist eine Konstante (unabhängig von und ) und kann durch Skalieren der Wahrscheinlichkeitsfunktion so gefunden werden, dass sie auf 1 integriert wird.$\frac{P(a,b)}{P(d)}$$a$$b$

Die Log-Likelihood-Funktion für iid-Variablen lautet:$n$$d_i \sim GN(0,a,b)$

Für diese Funktion sollte es nicht zu schwierig sein, sie zu zeichnen und ein Maximum zu finden.

Ich verstehe nicht ganz, wie Sie das Modell aufbauen: Insbesondere scheint es mir, dass für einen bestimmten Samen die möglichen Pollenausbreitungsentfernungen eine endliche Menge sind und Ihre "Ausbreitungswahrscheinlichkeit" daher besser als " Ausbreitungsrate "(da sie durch Summieren über mutmaßliche Väter normalisiert werden müsste, um eine Wahrscheinlichkeit zu sein). Daher haben die Parameter möglicherweise nicht ganz die Bedeutung (wie in plausiblen Werten), die Sie erwarten.

Ich habe in der Vergangenheit an einigen ähnlichen Problemen gearbeitet und werde daher versuchen, die Lücken in meinem Verständnis zu schließen, um einen möglichen Ansatz / einen kritischen Blick vorzuschlagen. Entschuldigung, wenn ich den Punkt Ihrer ursprünglichen Frage völlig verpasse. Die folgende Behandlung folgt im Wesentlichen Hadfield et al. (2006) , einem der besseren Artikel über diese Art von Modell.

Lassen des Genotyps an locus bezeichnen für einzelne . Für Nachkommen mit bekannter Mutter und mutmaßlichem Vater sei die Wahrscheinlichkeit der beobachteten Nachkommengenotypen - im einfachsten Fall ist dies einfach ein Produkt der Mendelschen Vererbungswahrscheinlichkeiten, aber in komplizierteren Fällen kann ein Modell für Genotypisierungsfehler oder fehlende elterliche Genotypen enthalten sein, daher füge ich Störparameter hinzu ) .$X_{l,k}$$l$$k$$i$$m_i$$f$

Sei der Pollenverteilungsabstand für die Nachkommen und sei der Abstand zwischen der bekannten Mutter und dem mutmaßlichen Vater und sei ist die Ausbreitungsrate (z. B. eine gewichtete Kombination von verallgemeinerten normalen und einheitlichen PDFs wie in Ihrer Frage). Um die Ausbreitungsrate als Wahrscheinlichkeit auszudrücken, normalisieren Sie wrt auf den endlichen Zustandsraum: die (endliche) Menge möglicher Ausbreitungsentfernungen, die durch die (endliche) Anzahl mutmaßlicher Väter in Ihrem Untersuchungsgebiet induziert werden, so dass $\delta_i$$i$$d_{m_i,f}$$m_i$$f$$D_{i,f} = q(d_{m_i,f}|a,b,c)$

Sei die väterliche Zuordnung von Samen , wenn Pflanze der Vater der Nachkommen . Unter der Annahme einer Uniform vor Vaterschaftszuweisungen ist Mit anderen Worten, abhängig von anderen Parametern und Genotypen ist die väterliche Zuordnung ein diskreter RV mit endlicher Unterstützung, der durch Integration über diese Unterstützung (mögliche Väter) normalisiert wird.$P_i$$i$$P_i = f$$f$$i$

Eine vernünftige Möglichkeit, einen einfachen Sampler für dieses Problem zu schreiben, ist Metropolis-in-Gibbs:

1. Unter der Bedingung, dass Vaterschaftszuweisungen für alle aktualisieren . Dies ist ein diskretes Wohnmobil mit endlicher Unterstützung, sodass Sie leicht eine exakte Probe zeichnen können$\{a,b,c,\theta\}$$P_i$$i$
2. Unter der Bedingung aktualisieren Sie mit einem Metropolis-Hastings-Update. Um das Ziel zu bilden, müssen nur die Werte in den obigen Gleichungen aktualisiert werden, sodass dies nicht kostspielig ist$\{P_i,\theta\}$$a,b,c$$D$
3. Conditional auf , update mit einem MH - Update. Um das Ziel zu bilden, müssen die Werte aktualisiert werden, was kostspielig ist, die Werte jedoch nicht.$\{P_i,a,b,c\}$$\theta$$G$$D$

Um die Kosten für das Zeichnen von Mustern von zu senken , können Sie die Schritte 1-2 vor 3 mehrmals ausführen. Um die Angebotsverteilungen in den Schritten 2-3 zu optimieren, können Sie Muster aus einem Vorlauf bis verwenden Schätzen Sie die Kovarianz der posterioren Gelenkverteilung für . Verwenden Sie dann diese Kovarianzschätzung in einem multivariaten Gaußschen Vorschlag. Ich bin sicher, dass dies nicht der effizienteste Ansatz ist, aber er ist einfach zu implementieren.$\{a,b,c\}$$\{a,b,c,\theta\}$

Nun, dieses Schema könnte nahe an dem liegen, was Sie bereits tun (ich kann anhand Ihrer Frage nicht sagen, wie Sie die Vaterschaft modellieren). Abgesehen von rechnerischen Bedenken ist mein größerer Punkt, dass die Parameter möglicherweise nicht die Bedeutung haben, von der Sie glauben, dass sie sie in Bezug auf die mittlere Ausbreitungsentfernung haben. Dies liegt daran , dass im Rahmen des Vaterschaftsmodell I oben beschrieben ist , zeigen sowohl treten in Zähler und Nenner (Normierungskonstante): so die räumliche Anordnung der Pflanzen einen potenziell starken Einfluss darauf haben, welche Werte von eine hohe Wahrscheinlichkeit oder eine hintere Wahrscheinlichkeit haben. Dies gilt insbesondere dann, wenn die räumliche Verteilung der Pflanzen ungleichmäßig ist.$a,b,c$$\mathrm{Pr}(P_i|\cdot)$$a,b,c$a , b , c$a,b,c$

Schließlich schlage ich vor, dass Sie sich das oben verlinkte Hadfield-Papier und das zugehörige R-Paket ("MasterBayes") ansehen, falls Sie dies noch nicht getan haben. Zumindest kann es Ideen liefern.