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:

1. Teraz przy użyciu tolerancji 0,001 dla średniej (było 1)
2. Zwiększono liczbę kroków do 500 tys., Aby uwzględnić mniejszą tolerancję
3. Zmniejszono SD rozkładu propozycji do 1, aby uwzględnić mniejszą tolerancję (było 10)
4. Dodano proste prawdopodobieństwo rnorm z n = 2k dla porównania
5. 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.

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$

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

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$

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.

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.