Продолжение про хиральность кошек

[dimview]

Объявленное ранее исследование продолжается.



Пока результат неопределённый:



Нужно больше данных.


Алгоритм:

• Распечатываем картинку
  
• Ставим точку в нижний левый угол, где координаты (0,0)
  
• Увидев кота, свернувшегося по часовой стрелке (вот так ⟳), смещаемся на единицу вправо
  
• Увидев кота, свернувшегося против часовой стрелки (вот так ⟲), смещаемся на единицу вверх
  
• Если добрались до красной линии, прекращаем тест. Кот предпочитает ⟲
  
• Если добрались до зелёной линии, прекращаем тест. Кот предпочитает ⟳
  
• Если добрались до жёлтой линии, прекращаем тест. Результат теста не является статистически значимым
  

Нулевая гипотеза — что коту всё равно и он сворачивается в обоих направлениях с одинаковой вероятностью p₀=0.5.

Альтернативная гипотеза — что коту не всё равно, то есть вероятность сворачиваться в одном направлении p₁ или меньше.

Мы хотим обнаружить эффект такого размера или больше (если он есть) с вероятностью 1-β.

Если эффекта нет, то мы хотим его найти с вероятностью меньше α. Оба направления возможны, поэтому у нас два хвоста (речь про хвосты распределения, а не про кошачьих хвосты), под каждым хвостом вероятность меньше α/2.

Проще всего выбрать фиксированное количество наблюдений (например, N=199) и критические значения. Выглядеть это будет так:



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

Начинаем в точке (0,0), смещаемся вправо на единицу если ⟳, смещаемся вверх на единицу если ⟲, продолжаем до тех пор, пока не упрёмся в одну из ограничительных линий, которые определяют результат. Синия тропинка показывает пример такого движения.

Тест с фиксированным количеством наблюдений требует ждать до конца. Если его остановить до того, как все N наблюдений собраны, то результат нельзя будет интерпретировать так, как было заранее задумано. Даже если мы видим 100 раз подряд ⟳ и ни одного ⟲, всё равно надо продолжать.

Поэтому будем использовать последовательный тест с ранней остановкой, похожий на тот, что описывает Evan Miller. У нас нет контрольной ячейки, поэтому использовать калькулятор оттуда не получится, но несложно посчитать отсечки самостоятельно и заодно помонтекарлить в R для проверки.

Нам нужно две отсечки, M и D. Для α=0.05, β=0.2, p₀=0.5 и p₁=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

Comments

[dimview]

Терминология заимствована у xaxam здесь.

Если бы у кошки было две стороны, как у монеты, а в других направлениях она была бы симметричной, тогда да. Но тут налицо весьма трёхмерный объект, у которого нос от хвоста отличается и лапы от спины тоже.

И вообще я тут наукой занимаюсь. Слово "хиральность", оно ж высоконаучное. А "на левом боку" — это ненаучный термин, потому что его все понимают.

[juan_gandhi]

В данном случае как раз хиральности нету. Можно взять, скажем, гвоздь, свернуть в колечко, и что? Вот оно лежит по солнышку; перевернули - лежит против солнышка. В три д нет проблем, если предмет можно перевернуть.

А вот с правилом буравчика уже не так.

Улыбнувшись

[dimview]

На фотографиях у кота нос и хвост находятся на разных уровнях (нос выше). Можно представить себе кота в виде одного витка резьбы, правой или левой. Гайка с правой резьбой на винт с левой резьбой не наворачивается, как её ни переворачивай.

Re: Улыбнувшись

[juan_gandhi]

Ха, так если кот повернет голову или хвост, у него хиральность меняется? Любопытно.

Утвердительно

[dimview]

Коты они такие, одним поворотом хвоста много чего могут.