Когда доверительный интервал «имеет смысл», а соответствующий достоверный интервал - нет?

Часто бывает так, что доверительный интервал с охватом 95% очень похож на вероятный интервал, который содержит 95% апостериорной плотности. Это происходит, когда предшествующий является однородным или почти однородным в последнем случае. Таким образом, доверительный интервал часто можно использовать для аппроксимации вероятного интервала и наоборот. Важно отметить, что из этого можно сделать вывод, что очень ошибочное неверное истолкование доверительного интервала как вероятного интервала практически не имеет практического значения для многих простых случаев использования.

Существует целый ряд примеров, когда этого не происходит, однако, похоже, что все они, похоже, подобраны сторонниками байесовской статистики в попытке доказать, что что-то не так с частым подходом. В этих примерах мы видим, что доверительный интервал содержит невозможные значения и т. Д., Что должно показывать, что они бессмысленны.

Я не хочу возвращаться к тем примерам или философскому обсуждению Байесовского против Фрикалистов.

Я просто ищу примеры обратного. Существуют ли случаи, когда доверительные и достоверные интервалы существенно различаются, а интервал, обеспечиваемый процедурой достоверности, явно превосходит?

Для пояснения: речь идет о ситуации, когда вероятный интервал обычно должен совпадать с соответствующим доверительным интервалом, то есть при использовании плоских, равномерных и т. Д. Априорных значений. Меня не интересует случай, когда кто-то выбирает произвольно плохой приор.

РЕДАКТИРОВАТЬ: В ответ на ответ @JaeHyeok Шин ниже, я должен не согласиться, что его пример использует правильную вероятность. Я использовал приблизительное байесовское вычисление для оценки правильного апостериорного распределения для тета, приведенного ниже в 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))

Это 95% вероятный интервал:

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

РЕДАКТИРОВАТЬ № 2:

Вот обновление после комментариев @JaeHyeok Shin. Я пытаюсь сделать это как можно более простым, но сценарий стал немного сложнее. Основные изменения:

Теперь используя допуск 0,001 для среднего (это было 1)
Увеличено количество шагов до 500 КБ, чтобы учесть меньший допуск
Уменьшил sd распределения предложения до 1, чтобы учесть меньший допуск (это было 10)
Добавлена простая вероятность rnorm с n = 2k для сравнения
Добавлен размер выборки (n) в качестве сводной статистики, установите допуск на 0,5 * n_target

Вот код:

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

Результаты, где hdi1 - это моя «вероятность», а hdi2 - простая норма (n, theta, 1):

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

Таким образом, после достаточного снижения допуска и за счет многих других шагов MCMC, мы можем увидеть ожидаемую ширину CrI для модели rnorm.

bayesian confidence-interval frequentist credible-interval синевато-багровый
источник

Не дублировать, но имеет тесную связь с stats.stackexchange.com/questions/419916/...

user158565

Как правило, когда у вас есть информативный априор, который является довольно неправильным, в неформальном смысле, например, Normal (0,1), когда фактическое значение равно -3,6, ваш вероятный интервал в отсутствие большого количества данных будет довольно плохим, когда смотрел с частой точки зрения.

jbowman

@jbowman Это особенно касается случая, когда используется равномерный априор или что-то вроде N (0, 1e6).

Livid

Несколько десятилетий назад настоящий байесовский статистик назвал статистику, которая ранее использовала неинформативный псевдо (или фальшивый) байесовский.

user158565

@ user158565 Это оффтоп, но единообразный априор - это всего лишь приближение. Если p (H_0) = p (H_1) = p (H_2) = ... = p (H_n), то все приоры могут выпасть из правил Байеса, что облегчает вычисления. Это не более неправильно, чем отбрасывать мелкие термины из знаменателя, когда это имеет смысл.

Livid

Ответы:

В последовательном экспериментальном проекте вероятный интервал может вводить в заблуждение.

(Отказ от ответственности: я не утверждаю, что это не разумно - это совершенно разумно в байесовских рассуждениях и не вводит в заблуждение с точки зрения байесовской точки зрения.)

$X$ $N(\theta,1)$ $\theta$ $n$ $\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\}.$

$P_{\theta}(N < \infty) = 1$ $\theta \in \mathbb{R}$

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

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

Возьмите домой сообщение: если вы заинтересованы в гарантии для частых, вы должны быть осторожны с использованием байесовских инструментов вывода, которые всегда действительны для байесовских гарантий, но не всегда для частых.

(Я узнал этот пример из удивительной лекции Ларри. Эта заметка содержит много интересных дискуссий о тонкой разнице между частыми и байесовскими системами. Http://www.stat.cmu.edu/~larry/=stat705/Lecture14.pdf )

РЕДАКТИРОВАТЬ В ABC Livid значение допуска слишком велико, поэтому даже для стандартной установки, в которой мы выбираем фиксированное количество наблюдений, оно не дает правильного CR. Я не знаком с ABC, но если я только изменил значение tol на 0,05, мы можем получить очень искаженный CR следующим образом

> 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

Конечно, цепь не очень хорошо стабилизирована, но даже если мы увеличим длину цепи, мы можем получить аналогичные CR - с перекосом в положительную часть.

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

JaeHyeok Shin
источник

«если мы установим достаточно большое значение k, то задняя часть θ будет приблизительно равна N (X_N, 1 / N)» . Мне кажется, что, очевидно, Pr (X | theta)! = Normal (theta, 1). То есть это неправильная вероятность для процесса, который породил вашу последовательность. Также есть опечатка. В исходном примере выборка прекращается, когда sqrt (n) * abs (mean (x))> k.

Livid

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

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

Пожалуйста, смотрите мои изменения в вопросе. Я все еще думаю, что ваш вероятный интервал не имеет смысла, потому что он использует неправильную вероятность. При использовании правильной вероятности, как в моем коде, мы получаем разумный интервал.

Livid

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$

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

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

Поскольку вероятный интервал формируется из апостериорного распределения на основе оговоренного предварительного распределения, вы можете легко построить очень плохой вероятный интервал, используя предварительное распределение, которое сильно сконцентрировано на крайне неправдоподобных значениях параметров. Вы можете создать надежный интервал, который «не имеет смысла», используя предыдущее распределение, полностью сконцентрированное на невозможных значениях параметров.

Восстановить Монику
источник

Или, что еще лучше, заслуживающий доверия, созданный настоятелем, который не согласен с вашим настоятелем (даже если это чужой настоятель) имеет хорошие шансы быть нелепым для вас. Это не редкость в науке; У меня были исследователи, которые говорили, что не хотят включать мнение экспертов, потому что в своих наблюдениях эксперты всегда были сильно самоуверенны.

Клифф А.Б.

Это конкретно о равномерных или "плоских" приорах.

Livid

@Livid: Вы должны обязательно включить, что вы говорите о плоских приорах в вашем вопросе. Это полностью меняет все.

Клифф AB

@CliffAB Это в первых двух предложениях, но я уточню, спасибо.

Livid

Если мы используем априорный флет, то это просто игра, в которой мы пытаемся придумать априорный флет при повторной калибровке, которая не имеет смысла.

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

Вот почему многие байесовцы возражают против плоских приоров.

Клифф AB
источник

Я объяснил свою мотивацию довольно четко. Я хочу что-то вроде примеров, где доверительные интервалы включают невозможные значения, но где вероятный интервал ведет себя хорошо. Если ваш пример основан на том, что вы делаете что-то бессмысленное, например, выбираете неправильную вероятность, то почему это кому-то интересно?

Livid

@Livid: функция правдоподобия вполне разумна. Плоским перед на лог-шансов нет. И это полный аргумент, который Байесов использует, чтобы сказать, что вы не должны использовать плоские приоры; на самом деле они могут быть чрезвычайно информативными и зачастую не такими, как задумал пользователь!

Клифф А.Б.

Вот Эндрю Гельман обсуждает некоторые вопросы плоских приоров.

Клифф AB

«Флэт-априор по лог-коэффициентам нет». Я имел в виду, что ставить фиксированную ставку на трансформированные логарифмические шансы для вас, похоже, бессмыслица, как ошибочное использование вероятности. Извините, но я не знаком с этим примером. Что именно должна делать эта модель?

Livid

@Livid: это может показаться необычным, но на самом деле это не так! Например, логистическая регрессия обычно учитывает все параметры в шкале логарифмов. Если бы вы использовали фиктивные переменные для всех ваших групп и использовали плоские априоры в параметрах регрессии, вы столкнулись бы именно с этой проблемой.

Клифф А.Б.