Rozkład propozycji uogólnionego rozkładu normalnego

Modeluję rozproszenie roślin przy użyciu uogólnionego rozkładu normalnego ( wpis na Wikipedii ), który ma funkcję gęstości prawdopodobieństwa:

gdzie jest przebytą odległością, jest parametrem skali, a jest parametrem kształtu. Średnią przebytą odległość podaje standardowe odchylenie tego rozkładu:$d$$a$$b$

Jest to wygodne, ponieważ pozwala na uzyskanie kształtu wykładniczego, gdy , kształtu Gaussa, gdy , i rozkładu lepeptycznego, gdy . Ta dystrybucja pojawia się regularnie w literaturze dotyczącej dyspersji roślin, chociaż ogólnie jest dość rzadka i dlatego trudno jest znaleźć informacje na jej temat.$b=1$$b=2$$b<1$

Najbardziej interesującymi parametrami są i średnia odległość rozproszenia.$b$

Próbuję oszacować i używając MCMC, ale walczę wymyślić skuteczny sposób, aby próbki wartości wniosku. Do tej pory korzystałem z Metropolis-Hastings i czerpałem z równomiernych rozkładów i , i otrzymuję średnie średnie odległości rozproszenia około 200-400 metrów, co ma sens biologiczny. Jednak konwergencja jest bardzo powolna i nie jestem przekonany, że bada pełną przestrzeń parametrów.$a$$b$$0 < a < 400$$0 < b<3$

Jego trudne wymyślić lepszej dystrybucji Wniosek dotyczący i , ponieważ zależą one od siebie, bez większego znaczenia na własną rękę. Średnia odległość rozproszenia ma wyraźne znaczenie biologiczne, ale daną średnią odległość rozproszenia można wyjaśnić nieskończenie wieloma kombinacjami i . Jako taki i są skorelowane w tylnej.$a$$b$$a$$b$$a$$b$

Do tej pory korzystałem z Metropolis Hastings, ale jestem otwarty na każdy inny algorytm, który by tu działał.

Pytanie: Czy ktoś może zaproponować bardziej skuteczny sposób rysowania wartości propozycji dla i ?$a$$b$

Edycja: Dodatkowe informacje o systemie: Badam populację roślin wzdłuż doliny. Celem jest określenie rozkładu odległości między pyłkami między roślinami dawcami a roślinami, które zapylają. Dane, które posiadam to:

1. Lokalizacja i DNA każdego możliwego dawcy pyłku
2. Nasiona zebrane z próbki 60 roślin matczynych (tj. Biorców pyłku), które zostały wyhodowane i genotypowane.
3. Lokalizacja i DNA każdej rośliny matczynej.

Nie znam tożsamości roślin dawcy, ale można to wywnioskować z danych genetycznych poprzez określenie, którzy dawcy są ojcami każdej sadzonki. Powiedzmy, że ta informacja jest zawarta w macierzy prawdopodobieństw G z wierszem dla każdego potomstwa i kolumną dla każdego kandydata na dawcę, co daje prawdopodobieństwo, że każdy kandydat będzie ojcem każdego potomstwa na podstawie wyłącznie danych genetycznych. Obliczenie G zajmuje około 3 sekund i musi być przeliczane przy każdej iteracji, co znacznie spowalnia.

Ponieważ ogólnie oczekujemy, że bliżsi kandydaci na dawców będą bardziej prawdopodobnie ojcami, wnioskowanie o ojcostwo jest dokładniejsze, jeśli wspólnie wnioskujesz o ojcostwie i rozproszeniu. Macierz D ma takie same wymiary jak G i zawiera prawdopodobieństwa ojcostwa oparte jedynie na funkcji odległości między matką a kandydatem i pewnym wektorze parametrów. Pomnożenie elementów w D i G daje wspólne prawdopodobieństwo ojcostwa na podstawie danych genetycznych i przestrzennych. Iloczyn pomnożonych wartości daje prawdopodobieństwo modelu rozproszenia.

Jak opisano powyżej, używam GND do modelowania rozproszenia. W rzeczywistości użyłem mieszanki GND i jednolitego rozkładu, aby umożliwić bardzo odległym kandydatom posiadanie większego prawdopodobieństwa ojcostwa z powodu samego przypadku (genetyka jest nieuporządkowana), który wzmógłby pozorny ogon GND, gdyby został zignorowany. Prawdopodobieństwo odległości rozproszenia wynosi:$d$

gdzie jest prawdopodobieństwem odległości rozproszenia od GND, N jest liczbą kandydatów, a ( ) określa, jaki wkład GND ma w rozproszenie.$\Pr(d|a,b)$$c$$0< c <1$

Istnieją zatem dwa dodatkowe czynniki, które zwiększają obciążenie obliczeniowe:

1. Odległość rozproszenia nie jest znana, ale należy ją wywnioskować przy każdej iteracji, a utworzenie G w tym celu jest kosztowne.
2. Istnieje trzeci parametr, , do zintegrowania.$c$

Z tych powodów wydawało mi się, że jest to trochę zbyt skomplikowane, aby wykonać interpolację siatki, ale cieszę się, że jestem przekonany, że jest inaczej.

Przykład

Oto uproszczony przykład kodu python, którego użyłem. Uprościłem szacowanie ojcostwa na podstawie danych genetycznych, ponieważ wymagałoby to wielu dodatkowych kodów i zastąpiłem je macierzą wartości od 0 do 1.

Najpierw zdefiniuj funkcje do obliczania 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

Następnie symuluj 2000 kandydatów i 800 potomstwa. Symuluj również listę odległości między matkami potomstwa i ojcami kandydującymi oraz obojętną macierzą G.

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]

Ustaw początkowe wartości parametrów:

# 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

Zaktualizuj kolejno a, b i c, a następnie oblicz współczynnik Metropolis.

# 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

Nie musisz używać metody Markov Chain Monte Carlo (MCMC).

Jeśli używasz wcześniejszych jednolitych rozkładów , robisz coś bardzo podobnego jak oszacowanie maksymalnego prawdopodobieństwa na ograniczonej przestrzeni dla parametrów i .$a$$b$

gdzie jest stałą (niezależną od i ) i można ją znaleźć poprzez skalowanie funkcji prawdopodobieństwa w taki sposób, że integruje się z 1.$\frac{P(a,b)}{P(d)}$$a$$b$

Funkcja wiarygodności dziennika dla zmiennych iid to:$n$$d_i \sim GN(0,a,b)$

Dla tej funkcji nie powinno być zbyt trudno wykreślić go i znaleźć maksimum.

Nie do końca rozumiem, w jaki sposób konfigurujesz model: w szczególności wydaje mi się, że dla danego materiału siewnego możliwe odległości rozproszenia pyłku są zbiorem skończonym, a zatem „prawdopodobieństwo rozproszenia” można lepiej określić jako „ współczynnik rozproszenia ”(ponieważ byłoby prawdopodobne, że trzeba go znormalizować poprzez zsumowanie przypuszczalnych ojców). Dlatego parametry mogą nie mieć takiego znaczenia (jak w prawdopodobnych wartościach), jakiego oczekujesz.

W przeszłości pracowałem nad kilkoma podobnymi problemami, dlatego spróbuję uzupełnić luki w moim zrozumieniu, aby zasugerować możliwe podejście / krytyczne spojrzenie. Przepraszam, jeśli całkowicie pomijam sedno twojego pierwotnego pytania. Poniższe leczenie jest zasadniczo zgodne z Hadfield i wsp. (2006) , jednym z lepszych artykułów na temat tego rodzaju modelu.

Niech oznacza genotyp w locus dla niektórych pojedynczych . Dla potomstwa ze znaną matką i domniemanym ojcem , niech prawdopodobieństwo zaobserwowanych genotypów potomstwa będzie - w najprostszym przypadku jest to po prostu iloczyn prawdopodobieństwa dziedziczenia przez Mendla, ale w bardziej skomplikowanych przypadkach może zawierać pewien model błędu genotypowania lub brakujących genotypów rodzicielskich, więc dołączam nieprzyjemny parametr (y) ) .$X_{l,k}$$l$$k$$i$$m_i$$f$

Niech będzie odległością rozproszenia pyłku dla potomstwa , i niech będzie odległością między znaną matką a domniemanym ojcem , i niech to współczynnik rozproszenia (np. ważona kombinacja uogólnionych normalnych i jednolitych plików pdf, jak w pytaniu). Aby wyrazić szybkość rozproszenia jako prawdopodobieństwo, unormuj wrt do skończonej przestrzeni stanów: (skończony) zestaw możliwych odległości rozproszenia indukowanych przez (skończoną) liczbę przypuszczalnych ojców w twoim obszarze badań, tak aby $\delta_i$$i$$d_{m_i,f}$$m_i$$f$$D_{i,f} = q(d_{m_i,f}|a,b,c)$

Niech będzie ojcowskim przypisaniem nasion , to znaczy jeśli roślina jest ojcem potomstwa . Zakładając jednolity wcześniejszy przydział na ojcostwo, .$P_i$$i$$P_i = f$$f$$i$

Tak więc rozsądnym sposobem napisania prostego samplera dla tego problemu jest Metropolis-within-Gibbs:

1. od zaktualizuj przypisania ojcostwa dla wszystkich . Jest to dyskretny rv ze skończonym wsparciem, dzięki czemu możesz łatwo narysować dokładną próbkę$\{a,b,c,\theta\}$$P_i$$i$
2. Warunkowo na zaktualizuj aktualizacją Metropolis-Hastings. Aby utworzyć cel, należy zaktualizować tylko wartości w powyższych równaniach, więc nie jest to kosztowne$\{P_i,\theta\}$$a,b,c$$D$
3. Warunkowo na , aktualizacja z aktualizacją MH. Aby utworzyć cel, wartości muszą zostać zaktualizowane, co jest kosztowne, ale nie.$\{P_i,a,b,c\}$$\theta$$G$$D$

Aby zmniejszyć koszt rysowania próbek , możesz wykonać kroki 1-2 wiele razy przed 3. Aby dostroić rozkłady propozycji w krokach 2-3, możesz użyć próbek ze wstępnego uruchomienia do oszacuj kowariancję wspólnego tylnego rozkładu dla . Następnie użyj tego oszacowania kowariancji w wielowymiarowej propozycji Gaussa. Jestem pewien, że nie jest to najbardziej wydajne podejście, ale jest łatwe do wdrożenia.$\{a,b,c\}$$\{a,b,c,\theta\}$

Teraz ten schemat może być zbliżony do tego, co już robisz (nie mogę powiedzieć, w jaki sposób modelujesz ojcostwo na podstawie twojego pytania). Ale poza kwestiami obliczeniowymi, moim większym punktem jest to, że parametry mogą nie mieć znaczenia, które według nich mają, w odniesieniu do średniej odległości rozproszenia. Wynika to z faktu, że w kontekście modelu ojcostwa który opisałem powyżej, wchodzą zarówno w licznik, jak i mianownik (stałą normalizującą): w ten sposób przestrzenny układ roślin będzie mieć potencjalnie silny wpływ na to, które wartości mają wysokie prawdopodobieństwo lub prawdopodobieństwo późniejsze. Jest to szczególnie prawdziwe, gdy rozkład przestrzenny roślin jest nierówny.$a,b,c$$\mathrm{Pr}(P_i|\cdot)$$a,b,c$a , b , c$a,b,c$

Na koniec, proponuję spojrzeć na ten dokument Hadfielda powiązany z powyższym i towarzyszący mu pakiet R („MasterBayes”), jeśli jeszcze tego nie zrobiłeś. Przynajmniej może dostarczyć pomysłów.