Kiedy przedział ufności „ma sens”, ale odpowiadający mu przedział wiarygodności nie?

Często zdarza się, że przedział ufności z 95% pokryciem jest bardzo podobny do wiarygodnego przedziału, który zawiera 95% gęstości tylnej. Dzieje się tak, gdy przeor jest jednolity lub prawie jednolity w tym drugim przypadku. Dlatego przedział ufności może być często stosowany do przybliżenia wiarygodnego przedziału i odwrotnie. Co ważne, możemy wywnioskować z tego, że bardzo zła interpretacja przedziału ufności jako wiarygodnego przedziału ma niewielkie lub żadne praktyczne znaczenie w wielu prostych przypadkach użycia.

Istnieje wiele przykładów przypadków, w których tak się nie dzieje, jednak wszystkie wydają się być wybrane przez zwolenników statystyk bayesowskich, próbując udowodnić, że coś jest nie tak z podejściem częstym. W tych przykładach widzimy, że przedział ufności zawiera wartości niemożliwe itp., Co ma wskazywać, że są nonsensowne.

Nie chcę wracać do tych przykładów ani filozoficznej dyskusji Bayesian vs. Frequentist.

Po prostu szukam przykładów przeciwnych. Czy istnieją przypadki, w których zaufanie i wiarygodne przedziały są zasadniczo różne, a przedział przewidziany w procedurze zaufania jest wyraźnie wyższy?

Wyjaśnienie: Chodzi o sytuację, w której zwykle oczekuje się, że wiarygodny przedział czasu zbiegnie się z odpowiednim przedziałem ufności, tj. Przy stosowaniu płaskich, jednolitych itp. Nie interesuje mnie sprawa, w której ktoś wybiera arbitralnie złego przeora.

EDYCJA: W odpowiedzi na poniższą odpowiedź @JaeHyeok Shina, nie mogę się zgodzić, że jego przykład wykorzystuje prawidłowe prawdopodobieństwo. Użyłem przybliżonego obliczenia bayesowskiego do oszacowania prawidłowego rozkładu tylnej dla theta poniżej w R:

### Methods ###
# Packages
require(HDInterval)

# Define the likelihood
like <- function(k = 1.2, theta = 0, n_print = 1e5){
  x    = NULL
  rule = FALSE
  while(!rule){
    x     = c(x, rnorm(1, theta, 1))
    n     = length(x)
    x_bar = mean(x)

    rule = sqrt(n)*abs(x_bar) > k

    if(n %% n_print == 0){ print(c(n, sqrt(n)*abs(x_bar))) }
  }
  return(x)
}

# Plot results
plot_res <- function(chain, i){
    par(mfrow = c(2, 1))
    plot(chain[1:i, 1], type = "l", ylab = "Theta", panel.first = grid())
    hist(chain[1:i, 1], breaks = 20, col = "Grey", main = "", xlab = "Theta")
}


### Generate target data ### 
set.seed(0123)
X = like(theta = 0)
m = mean(X)


### Get posterior estimate of theta via ABC ###
tol   = list(m = 1)
nBurn = 1e3
nStep = 1e4


# Initialize MCMC chain
chain           = as.data.frame(matrix(nrow = nStep, ncol = 2))
colnames(chain) = c("theta", "mean")
chain$theta[1]  = rnorm(1, 0, 10)

# Run ABC
for(i in 2:nStep){
  theta = rnorm(1, chain[i - 1, 1], 10)
  prop  = like(theta = theta)

  m_prop = mean(prop)


  if(abs(m_prop - m) < tol$m){
    chain[i,] = c(theta, m_prop)
  }else{
    chain[i, ] = chain[i - 1, ]
  }
  if(i %% 100 == 0){ 
    print(paste0(i, "/", nStep)) 
    plot_res(chain, i)
  }
}

# Remove burn-in
chain = chain[-(1:nBurn), ]

# Results
plot_res(chain, nrow(chain))
as.numeric(hdi(chain[, 1], credMass = 0.95))

Jest to 95% wiarygodny przedział:

> as.numeric(hdi(chain[, 1], credMass = 0.95))
[1] -1.400304  1.527371

EDYCJA 2:

Oto aktualizacja po komentarzach @JaeHyeok Shin. Staram się, aby było to tak proste, jak to możliwe, ale skrypt stał się nieco bardziej skomplikowany. Główne zmiany:

Teraz przy użyciu tolerancji 0,001 dla średniej (było 1)
Zwiększono liczbę kroków do 500 tys., Aby uwzględnić mniejszą tolerancję
Zmniejszono SD rozkładu propozycji do 1, aby uwzględnić mniejszą tolerancję (było 10)
Dodano proste prawdopodobieństwo rnorm z n = 2k dla porównania
Dodano wielkość próby (n) jako statystykę podsumowującą, ustaw tolerancję na 0,5 * n_cel

Oto kod:

### Methods ###
# Packages
require(HDInterval)

# Define the likelihood
like <- function(k = 1.3, theta = 0, n_print = 1e5, n_max = Inf){
  x    = NULL
  rule = FALSE
  while(!rule){
    x     = c(x, rnorm(1, theta, 1))
    n     = length(x)
    x_bar = mean(x)
    rule  = sqrt(n)*abs(x_bar) > k
    if(!rule){
     rule = ifelse(n > n_max, TRUE, FALSE)
    }

    if(n %% n_print == 0){ print(c(n, sqrt(n)*abs(x_bar))) }
  }
  return(x)
}


# Define the likelihood 2
like2 <- function(theta = 0, n){
  x = rnorm(n, theta, 1)
  return(x)
}



# Plot results
plot_res <- function(chain, chain2, i, main = ""){
    par(mfrow = c(2, 2))
    plot(chain[1:i, 1],  type = "l", ylab = "Theta", main = "Chain 1", panel.first = grid())
    hist(chain[1:i, 1],  breaks = 20, col = "Grey", main = main, xlab = "Theta")
    plot(chain2[1:i, 1], type = "l", ylab = "Theta", main = "Chain 2", panel.first = grid())
    hist(chain2[1:i, 1], breaks = 20, col = "Grey", main = main, xlab = "Theta")
}


### Generate target data ### 
set.seed(01234)
X    = like(theta = 0, n_print = 1e5, n_max = 1e15)
m    = mean(X)
n    = length(X)
main = c(paste0("target mean = ", round(m, 3)), paste0("target n = ", n))



### Get posterior estimate of theta via ABC ###
tol   = list(m = .001, n = .5*n)
nBurn = 1e3
nStep = 5e5

# Initialize MCMC chain
chain           = chain2 = as.data.frame(matrix(nrow = nStep, ncol = 2))
colnames(chain) = colnames(chain2) = c("theta", "mean")
chain$theta[1]  = chain2$theta[1]  = rnorm(1, 0, 1)

# Run ABC
for(i in 2:nStep){
  # Chain 1
  theta1 = rnorm(1, chain[i - 1, 1], 1)
  prop   = like(theta = theta1, n_max = n*(1 + tol$n))
  m_prop = mean(prop)
  n_prop = length(prop)
  if(abs(m_prop - m) < tol$m &&
     abs(n_prop - n) < tol$n){
    chain[i,] = c(theta1, m_prop)
  }else{
    chain[i, ] = chain[i - 1, ]
  }

  # Chain 2
  theta2  = rnorm(1, chain2[i - 1, 1], 1)
  prop2   = like2(theta = theta2, n = 2000)
  m_prop2 = mean(prop2)
  if(abs(m_prop2 - m) < tol$m){
    chain2[i,] = c(theta2, m_prop2)
  }else{
    chain2[i, ] = chain2[i - 1, ]
  }

  if(i %% 1e3 == 0){ 
    print(paste0(i, "/", nStep)) 
    plot_res(chain, chain2, i, main = main)
  }
}

# Remove burn-in
nBurn  = max(which(is.na(chain$mean) | is.na(chain2$mean)))
chain  = chain[ -(1:nBurn), ]
chain2 = chain2[-(1:nBurn), ]


# Results
plot_res(chain, chain2, nrow(chain), main = main)
hdi1 = as.numeric(hdi(chain[, 1],  credMass = 0.95))
hdi2 = as.numeric(hdi(chain2[, 1], credMass = 0.95))


2*1.96/sqrt(2e3)
diff(hdi1)
diff(hdi2)

Wyniki, w których hdi1 jest moim „prawdopodobieństwem”, a hdi2 to prosty rnorm (n, theta, 1):

> 2*1.96/sqrt(2e3)
[1] 0.08765386
> diff(hdi1)
[1] 1.087125
> diff(hdi2)
[1] 0.07499163

Po wystarczającym obniżeniu tolerancji i kosztem wielu kolejnych kroków MCMC, możemy zobaczyć oczekiwaną szerokość CrI dla modelu rnorm.

bayesian confidence-interval frequentist credible-interval Wściekły
źródło

Nie powiela się, ale ma ścisły związek ze stats.stackexchange.com/questions/419916/…

user158565

Ogólnie rzecz biorąc, jeśli masz uprzedzenie informacyjne, które jest całkiem błędne, w nieformalnym sensie, np. Normalny (0,1), gdy rzeczywista wartość wynosi -3,6, twój wiarygodny przedział przy braku dużej ilości danych będzie dość słaby, gdy patrząc z częstej perspektywy.

jbowman

@ jbowman Dotyczy to w szczególności przypadku używania munduru przed lub czegoś takiego jak N (0, 1e6).

Livid

Kilkadziesiąt lat temu prawdziwy Bayesian nazywał statystykiem, który używał nieinformacyjnego wcześniej jako pseudo- (lub fałszywy) Bayesian.

user158565,

@ user158565 To jest nie na temat, ale jednolity przeor jest tylko przybliżeniem. Jeśli p (H_0) = p (H_1) = p (H_2) = ... = p (H_n), wówczas wszystkie priorytety mogą zostać usunięte z reguły Bayesa, co ułatwia obliczenia. Nie jest to bardziej błędne niż usuwanie małych terminów z mianownika, gdy ma to sens.

Livid

Odpowiedzi:

W sekwencyjnym projekcie eksperymentalnym wiarygodny odstęp czasu może wprowadzać w błąd.

(Uwaga: nie twierdzę, że jest to nieuzasadnione - jest całkowicie uzasadnione w rozumowaniu bayesowskim i nie wprowadza w błąd z punktu widzenia bayesowskiego punktu widzenia.)

Dla prostego przykładu, powiedzmy, że mamy maszynę podającą nam losową próbkę $X$ z $N(\theta,1)$ o nieznanym $\theta$ . Zamiast rysunek $n$ próbek IId, zwracamy próbek aż $\sqrt{n} \bar{X}_n > k$ $k$ $N$

N = inf {n \geq 1 : \sqrt{n} {\bar{X}}_{n} > k} .

$N = \inf\left\{n \geq 1 : \sqrt{n}\bar{X}_n > k \right\}.$

Z prawa iterowanego logarytmu wiemy, że dla dowolnego . Ten typ reguły zatrzymywania jest powszechnie stosowany w sekwencyjnych testach / szacunkach w celu zmniejszenia liczby próbek do wnioskowania. $P_{\theta}(N < \infty) = 1$ $\theta \in \mathbb{R}$

Zasada prawdopodobieństwa pokazuje, że reguła zatrzymania nie wpływa na tylną stronę a zatem dla każdego rozsądnego gładkiego wcześniejszego (np. , jeśli ustawimy wystarczająco duże , tylna część ma w przybliżeniu a zatem wiarygodny przedział jest w przybliżeniu podawany jako Jednak z definicji wiemy, że ten wiarygodny przedział nie zawiera jeśli jest od tego czasu duże $\theta$ $\pi(\theta)$ $\theta \sim N(0, 10))$ $k$ $\theta$ $N(\bar{X}_N, 1/N)$

C I_{b a y e s} := [{\bar{X}}_{N} - \frac{1.96}{\sqrt{N}}, {\bar{X}}_{N} + \frac{1.96}{\sqrt{N}}] .

$CI_{bayes} :=\left[\bar{X}_N - \frac{1.96}{\sqrt{N}}, \bar{X}_N + \frac{1.96}{\sqrt{N}}\right].$

N

$N$

0

$0$

k

$k$

0 < {\bar{X}}_{N} - \frac{k}{\sqrt{N}} ≪ {\bar{X}}_{N} - \frac{1.96}{\sqrt{N}}

$0 < \bar{X}_N - \frac{k}{\sqrt{N}} \ll \bar{X}_N - \frac{1.96}{\sqrt{N}}$

dla . Dlatego częstość pokrycia wynosi zero, ponieważ a jest osiągane, gdy wynosi . Natomiast zasięg bayesowski jest zawsze równy w przybliżeniu ponieważ

k ≫ 0

$k \gg 0$

C I_{b a y e s}

$CI_{bayes}$

inf_{θ} P_{θ} (θ \in C I_{b a y e s}) = 0,

$\inf_{\theta} P_\theta( \theta \in CI_{bayes} ) = 0,$

0

$0$

θ

$\theta$

0

$0$

0.95

$0.95$

P (θ \in C I_{b a y e s} | X_{1}, \dots, X_{N}) \approx 0.95.

$\mathbb{P}( \theta \in CI_{bayes} | X_1, \dots, X_N) \approx 0.95.$

Weź wiadomość do domu: jeśli jesteś zainteresowany gwarancją częstego, powinieneś ostrożnie korzystać z narzędzi wnioskowania bayesowskiego, które zawsze obowiązują dla gwarancji bayesowskich, ale nie zawsze dla częstych.

(Nauczyłem się tego przykładu z niesamowitego wykładu Larry'ego. Ta notatka zawiera wiele interesujących dyskusji na temat subtelnej różnicy między frameworkami częstokrzyskimi a bayesowskimi. Http://www.stat.cmu.edu/~larry/=stat705/Lecture14.pdf )

EDYCJA W Livid's ABC wartość tolerancji jest zbyt duża, więc nawet dla standardowego ustawienia, w którym próbkujemy stałą liczbę obserwacji, nie daje prawidłowego CR. ABC nie jestem zaznajomiony z ABC, ale jeśli zmieniłem tylko wartość tol na 0,05, możemy mieć bardzo wypaczony CR jak poniżej

> X = like(theta = 0)
> m = mean(X)
> print(m)
[1] 0.02779672

> as.numeric(hdi(chain[, 1], credMass = 0.95)) [1] -0.01711265 0.14253673

Oczywiście łańcuch nie jest dobrze ustabilizowany, ale nawet jeśli zwiększymy długość łańcucha, możemy uzyskać podobny CR - przekrzywiony do części dodatniej.

W rzeczywistości myślę, że reguła odrzucania oparta na średniej różnicy nie jest dobrze dopasowana w tym ustawieniu, ponieważ z dużym prawdopodobieństwem jest bliski jeśli i bliski jeśli . $\sqrt{N}\bar{X}_N$ $k$ $0<\theta \ll k$ $-k$ $-k \ll \theta <0$

JaeHyeok Shin
źródło

„jeśli ustawimy wystarczająco dużą wartość k, tylna θ wynosi około N (X_N, 1 / N)” . Wydaje mi się, że oczywiście Pr (X | theta)! = Normal (theta, 1). To znaczy, że jest to niewłaściwe prawdopodobieństwo dla procesu, który wygenerował Twoją sekwencję. Jest też literówka. W oryginalnym przykładzie zatrzymujesz próbkowanie, gdy sqrt (n) * abs (mean (x))> k.

Livid,

Dziękuję za komentarze. Nawet pod regułą zatrzymania prawdopodobieństwo podaje się jako . Jest więc taki sam jak iloczyn N normalnych obserwacji, chociaż N jest teraz losowy. Ten przykład jest nadal aktualny dla obecnej reguły zatrzymania, chociaż ten, na który wskazałeś, jest oryginalnym i historycznym przykładem - Istnieje ciekawa historia tego, jak częsty i Bayesian debatowali z tą regułą zatrzymania. Możesz zobaczyć na przykład errorstatistics.com/2013/04/06/…

\prod_{i = 1}^{N} ϕ (X_{i} - θ)

$\prod_{i=1}^N \phi(X_i - \theta)$

JaeHyeok Shin

Proszę zobaczyć moją edycję w pytaniu. Nadal uważam, że Twój wiarygodny przedział czasu nie ma sensu, ponieważ wykorzystuje nieprawidłowe prawdopodobieństwo. Używając prawidłowego prawdopodobieństwa, jak w moim kodzie, otrzymujemy rozsądny odstęp czasu.

Livid,

Dzięki za szczegółowy eksperyment. W twoim ustawieniu jest zbyt małe, aby nierówności . Zauważ, że musi być większe niż 1,96, i dla porównania błędu aproksymacji, myślę, że byłoby bezpiecznym wyborem.

k

$k$

0 < {\bar{X}}_{N} - k / \sqrt{N} ≪ {\bar{X}}_{N} - 1.96 / \sqrt{N}

$0 < \bar{X}_N - k / \sqrt{N} \ll \bar{X}_N - 1.96 / \sqrt{N}$

k

$k$

k > 10

$k > 10$

JaeHyeok Shin

Czy możesz również ponownie sprawdzić, czy to obliczenie ABC działa w przypadku standardowym? Zastąpiłem funkcję \ code {like}, aby narysować stałą liczbę obserwacji (n = 2000). Następnie teoretyczna długość CR powinna wynosić około ale zawsze otrzymuję bardzo szerszy CR (długość wynosi około 2). Myślę, że poziom tolerancji jest zbyt duży.

2 \times 1.96 / \sqrt{2000} = 0.0876

$2 \times 1.96 / \sqrt{2000} = 0.0876$

JaeHyeok Shin

Ponieważ wiarygodny przedział jest tworzony z rozkładu tylnego, w oparciu o ustalony wcześniejszy rozkład, możesz łatwo skonstruować bardzo zły wiarygodny przedział, używając wcześniejszego rozkładu, który jest silnie skoncentrowany na wysoce nieprawdopodobnych wartościach parametrów. Można stworzyć wiarygodny przedział, który nie ma „sensu”, stosując wcześniejszy rozkład całkowicie skoncentrowany na niemożliwych wartościach parametrów.

Przywróć Monikę
źródło

Albo jeszcze lepiej: wiarygodna skonstruowana przez przeora, który nie zgadza się z twoim przeorem (chociaż jest to przeor kogoś innego), ma duże szanse na niedorzeczność. Nie jest to rzadkie w nauce; Powiedziałem naukowcom, że nie chcą uwzględniać opinii ekspertów, ponieważ w swoich obserwacjach eksperci zawsze byli zbyt pewni siebie.

Cliff AB,

Dotyczy to przede wszystkim jednolitych lub „płaskich” priorów.

Livid,

@Livid: Zdecydowanie powinieneś uwzględnić, że mówisz o płaskich priors w swoim pytaniu. To całkowicie zmienia wszystko.

Cliff AB,

@CliffAB Jest w pierwszych dwóch zdaniach, ale wyjaśnię, dzięki.

Livid

Jeśli używamy płaskiego przed, to jest to po prostu gra, w której staramy się wymyślić płaskie przed w parametryzacji, która nie ma sensu.

$\{0,1\}$ $\{1\}$

Dlatego wielu Bayesian sprzeciwia się płaskim przeorom.

Cliff AB
źródło

Wyjaśniłem moją motywację dość jasno. Chcę czegoś w rodzaju przykładów, w których przedziały ufności zawierają wartości niemożliwe, ale gdzie wiarygodny przedział zachowuje się dobrze. Jeśli twój przykład opiera się na robieniu czegoś bezsensownego, jak np. Wybranie niewłaściwego prawdopodobieństwa, to dlaczego miałoby to być interesujące dla kogokolwiek?

Livid

@Livid: funkcja wiarygodności jest całkowicie uzasadniona. Mieszkanie przed w dzienniku-kurs nie jest. I to jest cały argument, którego Bayesian używa, by powiedzieć, że nie powinieneś używać płaskich priorów; mogą być bardzo pouczające i często nie tak, jak zamierzał użytkownik!

Cliff AB

Oto Andrew Gelman omawiający niektóre kwestie dotyczące przeorów płaskich.

Cliff AB

„Płaskie wcześniejsze od log-odds nie jest”. Miałem na myśli to, że postawienie płaskiego wyniku na podstawie transformowanych logarytmów wydaje się dla ciebie nonsensem, jak użycie niewłaściwego prawdopodobieństwa. Przepraszam, ale nie znam tego przykładu. Co dokładnie ten model ma zrobić?

Livid

@Livid: może się to wydawać niezwykłe, ale tak naprawdę nie jest! Na przykład regresja logistyczna zazwyczaj uwzględnia wszystkie parametry w skali logarytmicznej. Gdybyś miał fikcyjne zmienne dla wszystkich grup i używał płaskich priorytetów w parametrach regresji, napotkałbyś dokładnie ten problem.

Cliff AB,