Mam dwie 2 godziny danych GPS z częstotliwością próbkowania 1 Hz (7200 pomiarów). Dane są podane w formie , gdzie jest niepewnością pomiaru.$(X, X_\sigma, Y, Y_\sigma, Z, Z_\sigma)$$N_\sigma$

Kiedy wezmę średnią ze wszystkich pomiarów (np. Średnią wartość Z tych dwóch godzin), jakie jest jej odchylenie standardowe? Mogę oczywiście obliczyć odchylenie standardowe od wartości Z, ale pomijam fakt, że znane są niepewności pomiaru ...

Edycja: wszystkie dane pochodzą z tej samej stacji, a wszystkie współrzędne są mierzone co sekundę. Ze względu na konstelacje satelitów itp. Każdy pomiar ma inną niepewność. Celem mojej analizy jest znalezienie przemieszczenia spowodowanego zdarzeniem zewnętrznym (trzęsieniem ziemi). Chciałbym wziąć średnią z 7200 pomiarów (2 godziny) przed trzęsieniem ziemi i inną średnią z 2 godzin po trzęsieniu ziemi, a następnie obliczyć wynikową różnicę (na przykład wysokość). Aby określić standardowe odchylenie tej różnicy, muszę znać odchylenie standardowe dwóch średnich.

Podejrzewam, że poprzednie odpowiedzi na to pytanie mogą nieco odbiegać od normy. Wydaje mi się, że to, o co tak naprawdę pyta oryginalny plakat, można by przeformułować jako „biorąc pod uwagę serię pomiarów wektorowych: z oraz kowariancja pomiaru :i=1,2,3,. . . ,7200Ci=( X 2 σ , i 0 0 0 Y

$$\vec{\theta}_{i} = \left( \begin{array}{c} X_{i} \\ Y_{i} \\ Z_{i} \end{array}\right)$$ $i=1, 2, 3,...,7200$ $$C_{i} = \left( \begin{array}{ccc} X_{\sigma,i}^{2} & 0 & 0 \\ 0 & Y_{\sigma,i}^{2} & 0 \\ 0 & 0 & Z_{\sigma,i}^{2} \end{array} \right)$$ jak miałbym poprawnie obliczyć średnią ważoną kowariancji dla tej serii pomiarów wektorowych, a następnie, jak poprawnie obliczyć jej odchylenie standardowe? ”Odpowiedź na to pytanie można znaleźć w wielu podręcznikach specjalizujących się w statystykach dla nauk fizycznych Jednym z przykładów, który szczególnie mi się podoba, jest Frederick James, „Metody statystyczne w fizyce eksperymentalnej” , drugie wydanie, World Scientific, 2006, sekcja 11.5.2, „Łączenie niezależnych szacunków”, str. 323–324. Kolejny bardzo dobry, ale bardziej tekstem wprowadzającym, który opisuje ważone wariancją średnie obliczanie wartości skalarnych (w przeciwieństwie do pełnych wielkości wektorów przedstawionych powyżej) jest Philip R. Bevington i D. Keith Robinson, „Redukcja danych i analiza błędów dla nauk fizycznych ”, 3. wydanie, McGraw-Hill, 2003, Sekcja 4.1.x, „Ważenie danych - niepewności Nonuniform”. Ponieważ w tym przypadku pytanie plakatu zawierało ukośną macierz kowariancji (tj. Wszystkie elementy poza przekątną są zerowe), problem można właściwie podzielić na trzy indywidualne (tj. X, Y, Z) problemy średniej ważonej skalarnie, więc analiza Bevington i Robinson ma tutaj równie dobre zastosowanie.

Ogólnie rzecz biorąc, odpowiadając na pytania stackexchange.com, zwykle nie uważam za użyteczne przepakowywanie długich pochodnych, które zostały już wcześniej przedstawione w wielu podręcznikach - jeśli chcesz naprawdę zrozumieć materiał i zrozumieć, dlaczego odpowiedzi wyglądają tak jak oni, to naprawdę powinieneś po prostu przeczytać wyjaśnienia, które zostały już opublikowane przez autorów podręcznika. Mając to na uwadze, przejdę bezpośrednio do powtórzenia odpowiedzi, które inni już udzielili. Frederick James, ustawiając , średnia ważona to: a kowariancja średniej ważonej wynosi:$N=7200$

$$\vec{\theta}_{mean} = \left( \sum_{i=1}^{N} C_{i}^{-1} \right)^{-1} \left( \sum_{i=1}^{N} C_{i}^{-1} \vec{\theta_{i}} \right)$$ $$C_{mean} = \left( \sum_{i=1}^{N} C_{i}^{-1} \right)^{-1}$$ Ta odpowiedź jest całkowicie ogólna i będzie ważna bez względu na forma , nawet dla macierzy kowariancji pomiaru niepiśmiennych.$C_{i}$

Ponieważ zdarza się, że kowariancje pomiarowe są w tym konkretnym przypadku ukośne, analizy Bevington i Robinson można również wykorzystać do obliczenia średnich ważonych wariancją dla poszczególnych , i . Forma odpowiedzi skalarnej jest podobna do formy odpowiedzi wektorowej: a wariancja to lub równoważnie, i podobnie dla$X_{i}$$Y_{i}$$Z_{i}$

$$X_{mean} = \frac{\sum_{i=1}^{N} \frac{X_{i}}{X_{\sigma,i}^{2}}}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}}$$ $$X_{\sigma,mean}^{2} = \frac{1}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}}$$ Ymean,Yσ,meanZmean,Zσ,mea $$X_{\sigma,mean} = \sqrt{\frac{1}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}}}$$ $Y_{mean}, Y_{\sigma, mean}$i . Krótki wpis na Wikipedii, który również pojawia się w tej samej odpowiedzi dla przypadku o wartości skalarnej, jest dostępny tutaj .$Z_{mean}, Z_{\sigma, mean}$

Należy to łatwo rozwiązać za pomocą wnioskowania bayesowskiego. Znasz właściwości pomiarowe poszczególnych punktów w odniesieniu do ich prawdziwej wartości i chcesz wnioskować o średniej populacji i SD, które wygenerowały prawdziwe wartości. To jest model hierarchiczny.

Ponowne odtworzenie problemu (podstawy Bayesa)

Zauważ, że podczas gdy ortodoksyjne statystyki dają ci jedną średnią, w ramach bayesowskich otrzymujesz rozkład wiarygodnych wartości średniej. Np. Obserwacje (1, 2, 3) z SD (2, 2, 3) mogły zostać wygenerowane na podstawie szacunku maksymalnego prawdopodobieństwa wynoszącego 2, ale również średnio 2,1 lub 1,8, choć nieco mniej prawdopodobne (biorąc pod uwagę dane) niż MLE. Więc oprócz SD, również wnioskujemy o średniej .

Kolejną różnicą pojęciową jest to, że przed dokonaniem obserwacji musisz zdefiniować swój stan wiedzy . Nazywamy to priors . Być może wiesz z góry, że określony obszar został zeskanowany i ma określony zakres wysokości. Całkowitym brakiem wiedzy byłoby posiadanie jednolitych (-90, 90) stopni jak wcześniej w X i Y i być może jednolitych (0, 10000) metrów wysokości (nad oceanem, poniżej najwyższego punktu na ziemi). Musisz zdefiniować rozkłady priors dla wszystkich parametrów, które chcesz oszacować, tj. Uzyskać rozkłady tylne . Dotyczy to również odchylenia standardowego.

Przeformułowując twój problem, zakładam, że chcesz wnioskować wiarygodne wartości dla trzech średnich (X.mean, Y.mean, X.mean) i trzech standardowych odchyleń (X.sd, Y.sd, X.sd), które mogłyby mieć wygenerował twoje dane.

Model

Używając standardowej składni BŁĘDÓW (do uruchomienia tego użyj pakietów WinBUGS, OpenBUGS, JAGS, stan lub innych), Twój model wyglądałby mniej więcej tak:

  model {
    # Set priors on population parameters
    X.mean ~ dunif(-90, 90)
    Y.mean ~ dunif(-90, 90)
    Z.mean ~ dunif(0, 10000)
    X.sd ~ dunif(0, 10)  # use something with better properties, i.e. Jeffreys prior.
    Y.sd ~ dunif(0, 10)
    Z.sd ~ dunif(0, 100)

    # Loop through data (or: set up plates)
    # assuming observed(x, sd(x), y, sd(y) z, sd(z)) = d[i, 1:6]
    for(i in 1:n.obs) {
      # The true value was generated from population parameters
      X[i] ~ dnorm(X.mean, X.sd^-2)  #^-2 converts from SD to precision
      Y[i] ~ dnorm(Y.mean, Y.sd^-2)
      Z[i] ~ dnorm(Z.mean, Z.sd^-2)

      # The observation was generated from the true value and a known measurement error
      d[i, 1] ~ dnorm(X[i], d[i, 2]^-2)  #^-2 converts from SD to precision
      d[i, 3] ~ dnorm(Y[i], d[i, 4]^-2)
      d[i, 5] ~ dnorm(Z[i], d[i, 6]^-2)
    }
  }

Oczywiście monitorujesz parametry .mean i .sd i używasz ich tylnych stron do wnioskowania.

Symulacja

Symulowałem niektóre takie dane:

# Simulate 500 data points
x = rnorm(500, -10, 5)  # mean -10, sd 5
y = rnorm(500, 20, 5)  # mean 20, sd 4
z = rnorm(500, 2000, 10)  # mean 2000, sd 10
d = cbind(x, 0.1, y, 0.1, z, 3)  # added constant measurement errors of 0.1 deg, 0.1 deg and 3 meters
n.obs = dim(d)[1]

Następnie uruchomiono model przy użyciu JAGS dla 2000 iteracji po wypaleniu 500 iteracji. Oto wynik dla X.sd.

Niebieski zakres wskazuje 95% przedział największej gęstości tylnej lub wiarygodnego (jeśli uważasz, że parametr jest obserwowany po zaobserwowaniu danych. Zauważ, że nie podaje tego zwykły przedział ufności).

Czerwona pionowa linia jest oszacowaniem MLE surowych danych. Zwykle jest tak, że najbardziej prawdopodobny parametr w estymacji Bayesa jest również najbardziej prawdopodobnym (maksymalnym prawdopodobieństwem) parametrem w statystykach ortodoksyjnych. Ale nie powinieneś zbytnio przejmować się górną częścią tylnej części ciała. Średnia lub mediana jest lepsza, jeśli chcesz sprowadzić ją do pojedynczej liczby.

Zauważ, że MLE / top nie ma wartości 5, ponieważ dane zostały wygenerowane losowo, a nie z powodu złych statystyk.

Ograniczenia

Jest to prosty model, który ma obecnie kilka wad.

1. Nie obsługuje tożsamości -90 i 90 stopni. Można to jednak zrobić, wprowadzając pewną zmienną pośrednią, która przesuwa skrajne wartości szacowanych parametrów do zakresu (-90, 90).
2. X, Y i Z są obecnie modelowane jako niezależne, chociaż prawdopodobnie są skorelowane i należy to wziąć pod uwagę, aby jak najlepiej wykorzystać dane. Zależy to od tego, czy urządzenie pomiarowe się poruszało (korelacja szeregowa i wspólny rozkład X, Y i Z dadzą wiele informacji), czy stoi w miejscu (niezależność jest w porządku). Na życzenie mogę rozwinąć odpowiedź, aby do tego podejść.

Powinienem wspomnieć, że istnieje wiele literatury na temat przestrzennych modeli bayesowskich, o których nie wiem.

Najpierw wprowadzam notację i ustawiam problem, korzystając z prostego podejścia, o którym wspomniałeś. Potem idź dalej. Użyję , aby zapoznać się z Z wektora dałeś.$\mathbf{z}$

Rozważ następujący model, w którym brakuje wyraźnego błędu pomiaru wzmianki: , gdzie jest szacunkową średnią wartością , a jest prawdziwą średnią wartością Z. Tutaj jest wektorem błędów w twoich danych i oczekujesz, że jeśli twoja próbka jest duża zbiegnie się w . Jeśli po prostu weźmiesz zaobserwowane wartości i uśrednisz je, otrzymasz a jeśli obliczysz standardowe odchylenie standardowe, otrzymasz , oszacowanie prawdziwego odchylenia standardowego populacjiˉ Z zμZϵ ˉ Z$\bar{Z} = \frac{\sum_{i=1}^n\mu_{Z} + \epsilon_i}{n}$$\bar{Z}$$\mathbf{z}$$\mu_Z$$\mathbf{\epsilon}$$\bar{Z}$$\mu_Z$$Z$$\bar{Z}$$\hat{\sigma}$$\sigma$ . Co jeśli chcesz skorzystać z wiedzy na temat błędu pomiaru?

Po pierwsze zauważmy, że możemy przeformułować początkowy model jako: , gdzie jest wektorem jedności, a skończy się na . Teraz to naprawdę wygląda jak regresja, ale nadal w zasadzie otrzymujemy oszacowanie . Jeśli wykonamy taką regresję, otrzymamy również oszacowanie standardowego błędu , który jest prawie tym, czego chcemy - to nic innego jak standardowy błąd (ale nadal chcemy uwzględnić błąd pomiaru).$\mathbf{z} = \mathbf{1}\beta + \epsilon$$\mathbf{1}$ˉ Z μ Z ϵ $\beta$$\bar{Z}$$\mu_Z$$\epsilon$$\mathbf{z}$

Możemy rozszerzyć nasz model początkowy, aby uzyskać model efektów mieszanych. , gdzie jest wektorem losowych efektów, a jest regresorem odnoszącym do . Podobnie jak w przypadku każdego efektu losowego, musisz przyjąć założenie dotyczące rozmieszczenia . Czy to prawda, że jest rozkładem błędu pomiaru dla$\mathbf{z} = \mathbf{1}\beta + \mathbf{Q}u +\epsilon$$u$oo U U Z σ $\mathbf{Q}$$\mathbf{z}$$u$$u$$Z_\sigma$$\mathbf{z}$? Jeśli tak, można to wykorzystać do zapewnienia rozkładu efektów losowych. Zazwyczaj oprogramowanie do wykonywania podstawowego modelowania efektów mieszanych zakłada, że ​​efekty losowe mają rozkład normalny (ze średnią 0 ...) i oszacuje dla ciebie wariancję. Być może możesz spróbować tego, aby przetestować koncepcję. Jeśli chcesz skorzystać z wcześniejszych informacji o rozkładzie błędu pomiaru, należy zastosować Bayesowski model efektów mieszanych. Możesz użyć R2OpenBUGS.

Po oszacowaniu tego modelu, standardowy błąd, który otrzymujesz dla resztek jest standardowym błędem, którym jesteś zainteresowany. Intuicyjnie, komponent efektów losowych modelu pochłania niektóre warianty, które możesz wyjaśnić, ponieważ wiesz, że istnieje pomiar błąd. Pozwala to uzyskać dokładniejsze oszacowanie wariantu$\epsilon$$\epsilon$

Zobacz ten artykuł, aby uzyskać głębszą dyskusję na temat tego podejścia losowych efektów w celu uwzględnienia błędu pomiaru. Twoja sytuacja jest podobna do tej, którą autorzy przedstawiają dla a jego błąd pomiaru jest zepsutej wersji . Przykład w sekcji 4 może dać wgląd w twoją sytuację.$\mathbf{D}$$\mathbf{W}$

Jak wspomniał Whuber, możesz chcieć uwzględnić autokorelację w swoich danych. Korzystanie z efektów losowych nie rozwiąże tego problemu.