О том как незнание центральной предельной теоремы позволяло отмывать бабло в средневековой Англии

Незнание того, как именно размер выборки влияет на статистические различия, создавало хаос на протяжении тысячелетия. В людском эпосе накопилось довольно много историй про это. Часть из них рассказана в статье “The Most Dangerous Equation”.

Что из себя может представлять опасное уравнение? Опасность может представлять два вида формул: те, которые вы знаете и те, которые вы не знаете. Первая категория уравнений может открывать для людей двери, за которыми находится смертельная опасность. Например, уравнение Эйнштейна открыло человечеству дорогу к атомной бомбе.

Однако нам интереснее другие уравнения. Уравнения, которые раскрывают свою опасность не тогда, когда мы знаем о них, а скорее тогда, когда мы их не знаем. Находясь под рукой, эти уравнения позволяют нам ясно понимать как работает природа. Незнание оставляет нас в опасном неведении.

В 1150 году было решено, что король Англии не может чеканить деньги и присваивать им любую ценность по своему выбору. Вместо этого ценность монеты должна была быть внутренней, основанной на том, сколько в её составе драгоценных материалов. Были установлены стандарты, сколько в какой монете должно быть золота. Чтобы проверить соответствует ли новая партия монет стандартам, проводили испытание Пикса.

Пикс (pyx) — это деревянный ящик, в котором находились стандартные монеты с правильным содержанием золота. Проверками занималась независимая организация, Монетный двор. Она состояла из баронов короля.

Монетный двор прекрасно понимал, что нельзя отчеканить абсолютно точную монету. Поэтому брали выборку новых монет, например сотню, и сравнивали её со стандартом. Надо было, чтобы эта сотня соответствовала заявленному уровню плюс минус некоторая погрешность. В качестве погрешности выбрали 1/400 веса.

Вес монет тогда измеряли в гранах, grain. Это единица измерения массы, основанная на весе ячменного зерна. Одна монета должна была весить 128 гран. Получается, что суммарный вес сотни монет должен был оказаться в диапазоне 12800 ± 32 гран.

Если мы пробили левую границу, надо казнить чеканщика за воровство. Если мы пробили правую границу, надо брать всю партию и переделывать, так как потрачено слишком много золота. Проблема в том, что интервал для суммы оказывался слишком широким, так как Монетный двор предполагал, что погрешность изменяется пропорционально числу монет. О том, что погрешность пропорциональна квадратному корню из числа монет, станет известно только через 600 лет благадоря Муавру.

Из-за слишком широкого интервала и страха казни, перекос был в сторону слишком тяжёлых монет. Скорее всего, были люди, которые чувствовали этот косяк в методологии Монетного двора и использовали его. Надо было просто брать из свежей партии монет самые тяжёлые, переплавлять их и навариваться на переплавке в типовые.

P.S. У автора оригинальной статьи почему-то 1/400 * 128 = 0.28. Видимо, он у баронов в доле.

Вот еще чем хотел поделиться

Вчера я закончил свой первый проект с использованием Typst (причем в Positron'е). До этого я просто экспериментировал с ним, а здесь вот небольшой проектик получился.

Не то чтобы в последнее время мне так часто нужен инструмент для верстки pdf-документов, но я никогда не испытывал особо теплых чувств к привычному многим LaTeX. А тут, новый инструмент, написанный на модном Rust. Да к тому же не просто язык разметки, а полноценный язык программирования (с переменными, функциями, управляющими конструкциями и т.д.). Я, честно сказать, в восторге, несмотря на достаточно сложные синтаксис и семантику. Полностью удовлетворен первым опытом.

Вообще, хочется сказать, что эти "хипстеры", пишущие на Rust, похоже, очень классные и полезные для общества люди (мечтаю к ним присоединиться, но уж очень сложный путь). Я периодически пробую использовать их альтернативы классическим инструментам, например, эмулятор терминала Rio, Zellij как замену Tmux, текстовый редактор Helix как замену Vim, GitUI для git и этот их великолепный менеджер проектов и пакетов Cargo, от всего только положительные эмоции.
Считаю, что им памятники из нержавеющей стали нужно ставить.

Для контекста можно посмотреть фильм "Человек, который познал бесконечность" и этот исторический анекдот от Райгородского

Кстати, Харди, который один из основных персонажей фильма и о котором Райгородский говорит, это вот тот самый Харди из формулы Харди-Вайнберга, которую, думаю, медики знают, помнят и любят

Блин, это гениально
https://www.youtube.com/watch?v=uqwC41RDPyg

Блин, вспомнил с этим Харди-Вайнбергом про свою медицинскую юность. Ностальгия прям. Полез посмотреть, что там сейчас с проектами, где я принимал участие.

Вот наткнулся, аж слезы на глазах )

2015 год, госпиталь Мулаго, Кампала, Уганда
Во втором ряду крайний слева – это я если что )

Здесь я этого почему-то не постил еще, исправляю.
Из курса МакЭлрита

И еще из закромов. Рубрика "Andrew Gelman says".

"...like Sander Greenland who disagrees with everybody"

В рекомендациях всплыло.

Если вы все еще думаете, что лекции по статистике – это скучно, то обязательно посмотрите эти выступления Kristin Lennox, она просто потрясающий рассказчик

Everything wrong with statistics (and how to fix it)
All About that Bayes: Probability, Statistics, and the Quest to Quantify Uncertainty

Качество не очень, но зато очень "залипательно"
https://www.youtube.com/watch?v=z23-awbMnGs

О, на блюскай прикольный проект затеяли

WeRateDAGs

Эту тему мы уже поднимали )

WeRateDAGs запостили

В этих постах, на мой взгляд, очень важная тема поднимается.
Там про сверхоптимизм при интерпретации размеров эффекта с очень необычным подходом к проблеме.

Там их 2, наверное, нужно последовательно читать, но я про второй как раз
https://www.linkedin.com/pulse/being-wrong-efficacy-probabilities-stephen-senn/?trackingId=p4y7N%2FJHQQyfLZM7YO9O9Q%3D%3D
https://www.linkedin.com/pulse/its-harder-than-you-think-stephen-senn/?trackingId=p4y7N%2FJHQQyfLZM7YO9O9Q%3D%3D

По следам сегодняшнего обсуждения в чате medstatistic решил сделать что-нибудь полезное.
Это будет байесовская лапласовская версия t-теста Уэлча*.

Почему я решил показать именно это (в смысле t-тест и именно такую его реализацию)? Просто потому что она мне самому кажется эстетически привлекательной. Представьте, люди очень давно могли посчитать на бумажке параметры апостериорного распределения (при строго определенных семействах априорных распределений, которые называются сопряженными с нормальным) для параметров нормального распределения матожидания (μ) и точности (1/σ²), но проблема в том, что это возможно сделать только если известны σ и μ, соответственно, то есть можно сказать это был абсолютно бесполезный навык с практической точки зрения.
Но вот в конце XX века появилась идея сэмплинга Гиббса – такого способа получения выборки из совместного распределения параметров, при котором мы на самом деле производим сэмплинг только из условных распределений. В данном случае, мы можем воспользоваться нашим древним знанием о сопряженных моделях для μ и 1/σ² и сделать очень простой, но при этом крайне эффективный сэмплер (демонстрация того, что это именно так приведена ниже).

* этим вариантом реализации не рекомендую пользоваться, ничего страшного не будет, но есть варианты лучше в смысле жесткости допущений о генеративном процессе

Итак, сэмплер Гиббса на основе сопряженных распределений:

normal_gibbs <- function(y, mu_0, tau_0, a, b, n_iter = 10^5) {
    n_obs <- length(y)
    sum_y <- sum(y)
    mu_sample <- tau_sample <- numeric(n_iter)
    
    # assigning starting values
    mu <- mean(y)
    tau <- 1 / var(y)
    
    # Gibbs sampling
    for (i in seq_len(n_iter)) {
        mu <- mu_sample[[i]] <- rnorm(
            n = 1,
            # parameters of conjugate posterior for μ given σ²
            mean = (sum_y * tau + mu_0 * tau_0) / (n_obs * tau + tau_0),
            sd = 1 / sqrt(n_obs * tau + tau_0)
        )
        tau <- tau_sample[[i]] <- rgamma(
            n = 1,
            # parameters of conjugate posterior for precision (1/σ²) given μ 
            shape = a + n_obs / 2,
            rate = b + 1 / 2 * sum((y - mu)^2)
        )
    }
    
    list("mu" = mu_sample, "sd" = sqrt(1/tau_sample))
}

Данные для примера:

set.seed(13799)

group_1 <- rnorm(10, mean = 120, sd = 15)
group_2 <- rnorm(15, mean = 130, sd = 25)

res_1 <- normal_gibbs(
    group_1,
    # parameters for prior for mu
    mu_0 = 100, tau_0 = 1 / (100^2),
    # parameters for prior for tau
    a = 0.01, b = 0.01
)

res_2 <- normal_gibbs(
    group_2,
    # parameters for prior for mu
    mu_0 = 100, tau_0 = 1 / (100^2),
    # parameters for prior for tau
    a = 0.01, b = 0.01
)