Продолжение про хиральность кошек
Jul. 28th, 2017 07:48 pm![[personal profile]](https://www.dreamwidth.org/img/silk/identity/user.png){kind=small}
dimviewОбъявленное ранее исследование продолжается.
Пока результат неопределённый:
Нужно больше данных.
Алгоритм:
Нулевая гипотеза — что коту всё равно и он сворачивается в обоих направлениях с одинаковой вероятностью p0=0.5.
Альтернативная гипотеза — что коту не всё равно, то есть вероятность сворачиваться в одном направлении p1 или меньше.
Мы хотим обнаружить эффект такого размера или больше (если он есть) с вероятностью 1-β.
Если эффекта нет, то мы хотим его найти с вероятностью меньше α. Оба направления возможны, поэтому у нас два хвоста (речь про хвосты распределения, а не про кошачьих хвосты), под каждым хвостом вероятность меньше α/2.
Проще всего выбрать фиксированное количество наблюдений (например, N=199) и критические значения. Выглядеть это будет так:
По горизонтали — сколько раз мы наблюдали свернувшегося по часовой стрелке кота с начала теста. По вертикали то же самое, только против часовой стрелки.
Начинаем в точке (0,0), смещаемся вправо на единицу если ⟳, смещаемся вверх на единицу если ⟲, продолжаем до тех пор, пока не упрёмся в одну из ограничительных линий, которые определяют результат. Синия тропинка показывает пример такого движения.
Тест с фиксированным количеством наблюдений требует ждать до конца. Если его остановить до того, как все N наблюдений собраны, то результат нельзя будет интерпретировать так, как было заранее задумано. Даже если мы видим 100 раз подряд ⟳ и ни одного ⟲, всё равно надо продолжать.
Поэтому будем использовать последовательный тест с ранней остановкой, похожий на тот, что описывает Evan Miller. У нас нет контрольной ячейки, поэтому использовать калькулятор оттуда не получится, но несложно посчитать отсечки самостоятельно и заодно помонтекарлить в R для проверки.
Нам нужно две отсечки, M и D. Для α=0.05, β=0.2, p0=0.5 и p1=0.4 получаем M=90 и D=33, соответственно картинка выглядит так:
Монте Карло намекает, что при наличии эффекта ранняя остановка уменьшает среднее число наблюдений до завершения теста:
Вот так можно подобрать значения M и D с помощью R.
Код строит два треугольника Паскаля, завёрнутые в ±D, один для нулевой гипотезы, другой для альтернативой, и увеличивает значения N и D до тех пор, пока не выполнятся требования к α и β.
В качестве примера вот треугольник для N=5 и D=3. Каждая строка — вектор вероятностей с индексом от 1 (который соответствует -D) до 2*D+1 (который соответствует +D). Первый и последний элементы вектора — вероятности того, что разница между cw и ccw достигнет ±D за N наблюдений.
Пока результат неопределённый:
Нужно больше данных.
Алгоритм:
- Распечатываем картинку
- Ставим точку в нижний левый угол, где координаты (0,0)
- Увидев кота, свернувшегося по часовой стрелке (вот так ⟳), смещаемся на единицу вправо
- Увидев кота, свернувшегося против часовой стрелки (вот так ⟲), смещаемся на единицу вверх
- Если добрались до красной линии, прекращаем тест. Кот предпочитает ⟲
- Если добрались до зелёной линии, прекращаем тест. Кот предпочитает ⟳
- Если добрались до жёлтой линии, прекращаем тест. Результат теста не является статистически значимым
Нулевая гипотеза — что коту всё равно и он сворачивается в обоих направлениях с одинаковой вероятностью p0=0.5.
Альтернативная гипотеза — что коту не всё равно, то есть вероятность сворачиваться в одном направлении p1 или меньше.
Мы хотим обнаружить эффект такого размера или больше (если он есть) с вероятностью 1-β.
Если эффекта нет, то мы хотим его найти с вероятностью меньше α. Оба направления возможны, поэтому у нас два хвоста (речь про хвосты распределения, а не про кошачьих хвосты), под каждым хвостом вероятность меньше α/2.
Проще всего выбрать фиксированное количество наблюдений (например, N=199) и критические значения. Выглядеть это будет так:
По горизонтали — сколько раз мы наблюдали свернувшегося по часовой стрелке кота с начала теста. По вертикали то же самое, только против часовой стрелки.
Начинаем в точке (0,0), смещаемся вправо на единицу если ⟳, смещаемся вверх на единицу если ⟲, продолжаем до тех пор, пока не упрёмся в одну из ограничительных линий, которые определяют результат. Синия тропинка показывает пример такого движения.
Тест с фиксированным количеством наблюдений требует ждать до конца. Если его остановить до того, как все N наблюдений собраны, то результат нельзя будет интерпретировать так, как было заранее задумано. Даже если мы видим 100 раз подряд ⟳ и ни одного ⟲, всё равно надо продолжать.
Поэтому будем использовать последовательный тест с ранней остановкой, похожий на тот, что описывает Evan Miller. У нас нет контрольной ячейки, поэтому использовать калькулятор оттуда не получится, но несложно посчитать отсечки самостоятельно и заодно помонтекарлить в R для проверки.
Нам нужно две отсечки, M и D. Для α=0.05, β=0.2, p0=0.5 и p1=0.4 получаем M=90 и D=33, соответственно картинка выглядит так:
- Держим два счётчика, cw и ccw, изначально оба равны нулю
- Увидев ⟳, увеличиваем cw на единицу
- Увидев ⟲, увеличиваем ccw на единицу
- Если cw - ccw == D, прекращаем тест. Кот предпочитает ⟳
- Если ccw - cw == D, прекращаем тест. Кот предпочитает ⟲
- Если cw == M или ccw == M, прекращаем тест. Результат не является статистически значимым
Монте Карло намекает, что при наличии эффекта ранняя остановка уменьшает среднее число наблюдений до завершения теста:
| Нет эффекта | Минимальный эффект | Худший случай | |
| Фиксированный размер | 199 | 199 | 199 |
| Ранняя остановка | 188 | 153 | 211 |
Вот так можно подобрать значения M и D с помощью R.
alpha <- 0.05 # Probability of finding an effect that does not exist
beta <- 0.2 # Probability of not finding an effect that exists
p0 <- 0.5 # Proportion under null hypothesis
p1 <- 0.4 # Proportion under alternative hypothesis (with minimal detectable effect)
give_up <- 10000 # Large number to fail early
for (d in 1:give_up) {
b0 <- numeric(2 * d + 1) # Vector of probabilities under null hypothesis
b0[d + 1] <- 1
b1 <- numeric(2 * d + 1) # Vector of probabilities under alternative hypothesis
b1[d + 1] <- 1
for (n in 1:give_up) { # n is number of observations
b0 <- c(b0[1], 0, b0[2:(length(b0)-1)]*p0) + c(b0[2:(length(b0)-1)]*(1-p0), 0, b0[length(b0)])
b1 <- c(b1[1], 0, b1[2:(length(b1)-1)]*p1) + c(b1[2:(length(b1)-1)]*(1-p1), 0, b1[length(b1)])
if ((b0[1] > alpha / 2) || (b1[1] > 1 - beta)) {
break;
}
}
if ((b0[1] < alpha / 2) && (b1[1] > 1 - beta)) {
break;
}
}
m <- (n - d) / 2 + 1
sprintf("m=%d, d=%d, fp=%f (should be <%f), tp=%f (should be >%f)",
m, d, b0[1], alpha / 2, b1[1], 1 - beta)
Код строит два треугольника Паскаля, завёрнутые в ±D, один для нулевой гипотезы, другой для альтернативой, и увеличивает значения N и D до тех пор, пока не выполнятся требования к α и β.
В качестве примера вот треугольник для N=5 и D=3. Каждая строка — вектор вероятностей с индексом от 1 (который соответствует -D) до 2*D+1 (который соответствует +D). Первый и последний элементы вектора — вероятности того, что разница между cw и ccw достигнет ±D за N наблюдений.
| N | 1 | 2 | 3 | 4 | 5 | 6 | 7 | ||||||
| 0 | 1 | ||||||||||||
| ↙ | ↘ | ||||||||||||
| 1 | .5 | .5 | |||||||||||
| ↙ | ↘ | ↙ | ↘ | ||||||||||
| 2 | .25 | .5 | .25 | ||||||||||
| ↙ | ↘ | ↙ | ↘ | ↙ | ↘ | ||||||||
| 3 | .125 | .375 | .375 | .125 | |||||||||
| ↓ | ↙ | ↘ | ↙ | ↘ | ↓ | ||||||||
| 4 | .125 | .1875 | .375 | .1875 | .125 | ||||||||
| ↓ | ↙ | ↘ | ↙ | ↘ | ↙ | ↘ | ↓ | ||||||
| 5 | .21875 | .28125 | .28125 | .21875 |
no subject
Date: 2017-07-28 06:45 pm (UTC)no subject
Date: 2017-07-28 07:41 pm (UTC)Если бы у кошки было две стороны, как у монеты, а в других направлениях она была бы симметричной, тогда да. Но тут налицо весьма трёхмерный объект, у которого нос от хвоста отличается и лапы от спины тоже.
И вообще я тут наукой занимаюсь. Слово "хиральность", оно ж высоконаучное. А "на левом боку" — это ненаучный термин, потому что его все понимают.
no subject
Date: 2017-07-28 07:59 pm (UTC)А вот с правилом буравчика уже не так.
Улыбнувшись
Date: 2017-07-28 08:37 pm (UTC)Re: Улыбнувшись
Date: 2017-07-28 09:05 pm (UTC)Утвердительно
Date: 2017-07-28 10:45 pm (UTC)no subject
Date: 2017-09-04 01:01 am (UTC)Спасибо за ссылку
Date: 2017-09-04 02:13 am (UTC)Слегка пугает то, что вышеупомянутая статья заканчивается фразой "those having a more intimate acquaintance with cattle will find it easy to give an answer." Вдруг за прошедшее с 1927 года время нашлись какие-нибудь британские учёные, готовые (разумеется, исключительно в научных целях) выполнить заветы авторов насчёт "более интимного знакомства со скотом". Страшно даже гуглить.