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:

\frac{b}{2 a Γ (1 / b)} e^{- (\frac{d}{a})^{b}}

$\frac{b}{2a\Gamma(1/b)} e^{-(\frac{d}{a})^b}$

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$

\sqrt{\frac{a^{2} Γ (3 / b)}{Γ (1 / b)}}

$\sqrt{\frac{a^2 \Gamma(3/b)}{\Gamma(1/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:

Lokalizacja i DNA każdego możliwego dawcy pyłku
Nasiona zebrane z próbki 60 roślin matczynych (tj. Biorców pyłku), które zostały wyhodowane i genotypowane.
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$

c Pr (d | a, b) + \frac{(1 - c)}{N}

$c \Pr(d|a,b) + \frac{(1-c)}{N}$

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:

Odległość rozproszenia nie jest znana, ale należy ją wywnioskować przy każdej iteracji, a utworzenie G w tym celu jest kosztowne.
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

distributions bayesian mcmc tellis
źródło

Szukasz a priori na A i B lub dystrybucji propozycji w algorytmie Metropolis-Hastings? Wygląda na to, że używałeś obu terminów zamiennie.

Robin Ryder,

Masz rację - przepraszam, że nie wyjaśniłeś. Najbardziej interesuje mnie dystrybucja propozycji dla MH. Zmieniłem tytuł, w którym odpowiednio wspomniałem o priors.

tellis

Pod mieszkaniem lub Jeffreys przed , tj. lub Uważam, że zmiana zmiennej na powoduje zamknięcie -form warunkowe .

a

$a$

π (a) \propto 1

$\pi(a)\propto 1$

π (a) \propto 1 / a

$\pi(a)\propto 1/a$

α = a^{- b}

$\alpha=a^{-b}$

π (a | b, data)

$\pi(a|b,\text{data})$

Xi'an

Zobacz en.wikipedia.org/wiki/Generalized_normal_distribution

kjetil b halvorsen

Pozostaje dość niejasne, czy jesteś zainteresowany ustaleniem przełożonego, który pomaga, czy wydajniejszym prowadzeniem Metropolis-Hastings.

Xi'an

Odpowiedzi:

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$

P (a, b; d) = P (d; a, b) \frac{P (a, b)}{P (d)} = L (a, b; d) \times c o n s t

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

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

\log L (a, b; d) = - n \log (2 a) - n \log (\frac{Γ (1 / b)}{b}) - \frac{1}{a^{b}} \sum_{i = 1}^{n} {(d_{i})}^{b}

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

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

Sextus Empiricus
źródło

Ta interpolacja wiązki działałaby dla dwóch parametrów i obserwowanych odległości i może być tym, co ostatecznie skończę. W rzeczywistości dokonuję wspólnej oceny odległości rozproszenia i wnioskowania o ojcostwie, która obejmuje co najmniej jeden inny parametr do zintegrowania, oraz dodatkowy okres prawdopodobieństwa, który jest naprawdę wolny (~ 3 sekundy na iterację), co naprawdę spowalnia łańcuch. Myślę, że potrzebowałbym około 10-krotnie więcej iteracji niż obecnie w przypadku łańcucha markowa.

tellis

@tellis te terminy „dystans rozproszenia” i „wnioskowanie o ojcostwie” nie do końca rozumiem. Może możesz podać bardziej konkretne informacje, dodając zestaw danych lub jego część. Robiąc to, możesz równie dobrze mówić o „jednym innym parametrze”. Więc, co dane jest to, że ty nie masz?

Sextus Empiricus

Dodałem przykład wykorzystujący dane symulowane.

tellis

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$

G_{i, f} = \prod_{l} P r (X_{l, i} | X_{l, m_{i}}, X_{l, f}, θ)

$G_{i,f} = \prod_l \mathrm{Pr}(X_{l,i}|X_{l,m_i}, X_{l,f}, \theta)$

θ

$\theta$

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

{\tilde{D}}_{i, f} = P r (δ_{i} = d_{m_{i}, f} | a, b, c) = \frac{D_{i, f}}{\sum_{k} D_{i, k}}

$\tilde{D}_{i,f} = \mathrm{Pr}(\delta_i = d_{m_i,f}|a,b,c) = \frac{D_{i,f}}{\sum_k D_{i,k}}$

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$

P r (P_{i} = f | a, b, c, θ, 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}}

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

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

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

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.

Nate Pope
źródło

Moje podejście jest rzeczywiście wzorowane na Hadfieldu, z dwiema istotnymi zmianami: (1) nasiona matki mogą być pełnym rodzeństwem, a zatem nie być niezależne. Problemem jest zatem wspólne wnioskowanie o strukturze rodzeństwa, ojcostwa, nietoperza. (2) Stosuję ułamkowe podejście do ojcostwa, aby rozważać wszystkich kandydatów jednocześnie proporcjonalnie do ich prawdopodobieństwa ojcostwa, zamiast aktualizować kolejno zadania ojcostwa, ponieważ istnieje duża przestrzeń do zbadania przez potencjalnych ojców.

tellis

Używam pakietu FAPS do robienia tych rzeczy.

tellis

Moje pytanie dotyczy w zasadzie wydajnego podziału propozycji na wykonanie punktu 2. Reszta twojej odpowiedzi opisuje coś bardzo zbliżonego do tego, co już zrobiłem, w tym sformułowanie produktu G i D (ale dzięki za to - nie byłem na pewno zrobiłem to poprawnie, więc warto wiedzieć, że druga para oczu się zgadza!).

tellis

Przykro mi, ale nie mam rozwiązania w zakresie dystrybucji propozycji w puszkach. Ale mam kilka spostrzeżeń: (1) Kroki 1-2 są bardzo tanie i można je powtarzać wiele razy niewielkim kosztem przed przejściem do kroku 3. Nawet z tandetną propozycją w kroku 2, wiele iteracji powinno wystarczyć do „ wykonywać duże ruchy ”w przestrzeni stanu .

a, b, c

$a,b,c$

Nate Pope

(2) Rozkład warunkowy w kroku 2 jest trójwymiarowy. Jak w: łatwy do wizualizacji. Jak wygląda nienormalizowany cel w oszacowaniu MAP przypisań ojcostwa dla ustalonego ? Wizualizacja nienormalizowanego celu dla różnych ojców powinna dać ci poczucie, czy jest multimodalny, płaski na obszarach itp.

a, b, c

$a,b,c$

G

$G$

Nate Pope