Применимость статистических тестов к задаче про саботажника

[realeugene]

Geen в сообщении #1645072 писал(а):

Я запустил тест 1000 раз и два раза выпало 2 нуля...

Тривиально считается. Вероятность того, что из 1000 тестов не менее двух раз выпадет ровно два нуля, составляет 0.09. Для 1000 тестов количество тестов с двумя нулями - пуассон с параметром $\mu=0.5$.

[Geen]

EUgeneUS в сообщении #1645194 писал(а):

А пока тест отверг нулевую гипотезу в 2-х случаях из 1000.

А ведь это с поправкой на множественный тест

[EUgeneUS]

Geen в сообщении #1645196 писал(а):

А ведь это с поправкой на множественный тест

Так у остальных 998-ми Вы получили p-value ровно единицу
Посмотрим, что будет без поправки

[Geen]

EUgeneUS в сообщении #1645197 писал(а):

получили p-value ровно единицу

Не p-value, а q-value...

[Geen]

EUgeneUS в сообщении #1645194 писал(а):

Посмотрим, что будет с поправкой.

Исправил по Вашему замечанию.

Используется синтаксис Python

def freq(arr,n):
        res = [0 for _ in range(n+1)];
        for _ in arr:   res[_] += 1;
        return res;
def bonferroni(arr):
        arr.sort();
        n = len(arr);
        res = reversed(list(_*n/(i+1) for i,_ in enumerate(arr)));
        cur = inf;
        res = [cur:=min(_,cur) for _ in res];
        res.reverse();
        return res;
p = 1/3;
n = 18;
m = 47;
k = 1000;
nn = floor(log10(n-1))+1;
kn = floor(log10(k-1))+1;
res = [[] for _ in range(n+1)];
pval_true = [];
pval_euus = [];
bin = [nchoosek(n,_)*pow(p,_)*pow(1-p,n-_)*m for _ in range(n+1)];
print(bin);
for _ in range(k):
        sample = list(sum(random()<p for i in range(n)) for i in range(m));
        frq = freq(sample,n);
        for i,_ in enumerate(frq):
                if _:
                        res[i].append(_);
        chi2_true = sum((frq[i]-bin[i])**2/bin[i] for i in range(n+1));
        pval_true.append(gammaincreg(n,chi2_true));
        chi2_euus = 0;
        cbin = bin.copy();
        dir = 0;
        tn = 0;
        while frq:
                cf = 0; cb = 0;
                while cf < 5 and frq:
                        cf += frq.pop(dir);
                        cb += cbin.pop(dir);
                chi2_euus += (cf-cb)**2/cb;
                tn += 1;
                dir = -1-dir;
        pval_euus.append(gammaincreg(tn-1,chi2_euus));
ma = max(max(_ or [0]) for _ in res);
for i,_ in enumerate(res):
        frq = freq(_,ma);
        #frq.pop(0);
        print(('%%%ii'%nn)%i,*[('%%%ii'%kn)%_ if _ else ' '*kn for _ in frq]);
qval_true = bonferroni(pval_true);
qval_euus = bonferroni(pval_euus);
print(qval_true[:10]);
print(qval_euus[:10]);
print(pval_true[:10]);
print(pval_euus[:10]);
 

Один из типичных результатов (как и ожидалось, стало существенно хуже "для теста"):

Используется синтаксис Text

[0.03180205577614669, 0.2862185019853202, 1.2164286334376107, 3.2438096891669606, 6.08214316718805, 8.515000434063268, 9.224583803568539, 7.906786117344461, 5.4359154556743166, 3.0199530309301754, 1.3589788639185787, 0.49417413233402846, 0.1441341219307583, 0.033261720445559605, 0.005939592936707071, 0.0007919457248942759, 7.424491170883835e-05, 4.367347747578727e-06, 1.2131521521052015e-07]
 0      30   2
 1     219  32   1
 2     377 212  93  37   4   1
 3     134 209 249 168 113  51  27   7   3   1
 4       9  35  79 130 141 190 132 109  86  37  30  15   4
 5       1   2  16  31  60  94 139 154 135 125  99  65  41  17  13   5   2   1
 6       1   2   8  23  63  74 119 142 141 140 117  62  48  27  18   9   5   1
 7       1   4  30  50  78 130 133 158 131 135  72  35  21  15   4   2       1
 8      24  52 117 188 161 176 123  69  42  27   8   5   1
 9     147 248 193 191 108  44  18   6   4   1   1
10     338 231 108  46   7   2
11     294  85  11   1
12     118  12   1
13      34
14       2
15       3
16
17
18
[0.0, 0.0, 0.0, 1.647080276493533e-50, 4.384365466863698e-50, 8.814699417233117e-35, 1.5946807957175212e-34, 9.489475524470615e-12, 4.087534256661161e-09, 5.995613895415072e-07]
[6.890291968253561e-05, 0.000634862728445439, 0.0054953010528217515, 0.0054953010528217515, 0.013633385613349416, 0.0281260422386174, 0.0281260422386174, 0.05944122196355882, 0.05944122196355882, 0.05976497793356464]    
[0.0, 0.0, 0.0, 6.588321105974132e-53, 2.192182733431849e-52, 5.28881965033987e-37, 1.1162765570022649e-36, 7.591580419576493e-14, 3.678780830995045e-11, 5.995613895415072e-09]
[6.890291968253562e-08, 1.269725456890878e-06, 2.112549143356511e-05, 2.1981204211287005e-05, 6.816692806674708e-05, 0.00018028249553454483, 0.0001968822956703218, 0.0005084917178444854, 0.0005349709976720294, 0.0007535425916394196]
 

-- 04.07.2024, 23:14 --

EUgeneUS
Вы ранее упоминали про "тяжёлые хвосты". Так вот, это хвосты не тяжёлые, а просто их распределение существенно отличается от нормального - поэтому Ваш тест на них и "врёт"...
Для таких перекошенных распределений 5 измерений на "бин" вовсе не достаточно...
И в частности поэтому, полсотни измерений совсем мало.

[EUgeneUS]

Geen
Еще просьба.
Можете посчитать, сколько раз из 1000 отвергается нулевая гипотеза по условию $p_value < 0.05$ при тесте с группировкой?

-- 04.07.2024, 23:21 --

Geen в сообщении #1645201 писал(а):

Так вот, это хвосты не тяжёлые, а просто их распределение существенно отличается от нормального - поэтому Ваш тест на них и "врёт"

ОМГ. Нам не важен вид ожидаемого распределения.

Завтра, если будет время, сделаю модель на 1000 кейсов в Excel. Посмотрим, насколько будет различаться с Вашим результатом.

[Geen]

EUgeneUS в сообщении #1645203 писал(а):

Можете посчитать, сколько раз из 1000 отвергается нулевая гипотеза по условию $p_value < 0.05$ при тесте с группировкой?

Это случайная величина с большой дисперсией (и распределением, которое я совершенно не знаю как посчитать).... Смысл?

[EUgeneUS]

FGJ, минимум 5 ожидаемых отсчетов на бин - это минимальные рекомендации, которые встречал. Но и наиболее частые.
Еще встречал такие:
1. Минимум 5 отсчетов на бин и минимум 100 всего.
2. Минимум 10 остчетов на бин, но для анализа таблиц совместимости $2 \times 2$
3. Минимум 10 отсчетов на бин, но может быть уменьшено до 5 или даже до 3 при количестве степеней свободы в несколько десятков.

-- 04.07.2024, 23:31 --

Geen в сообщении #1645204 писал(а):

Смысл?

Сравнить с ожидаемым значением $50$.

[Geen]

EUgeneUS в сообщении #1645205 писал(а):

... рекомендации ...

....
Если распределение в каждом бине близко к нормальному, то такую рекомендацию можно принять. Но в данном случае это не так.

EUgeneUS в сообщении #1645203 писал(а):

Можете посчитать, сколько раз из 1000 отвергается нулевая гипотеза по условию $p_value < 0.05$ при тесте с группировкой

Код

Используется синтаксис Python

        def freq(arr,n):
                res = [0 for _ in range(n+1)];
                for _ in arr:   res[_] += 1;
                return res;
        def bonferroni(arr):
                arr.sort();
                n = len(arr);
                res = reversed(list(_*n/(i+1) for i,_ in enumerate(arr)));
                cur = inf;
                res = [cur:=min(_,cur) for _ in res];
                res.reverse();
                return res;
        bigK = 100;
        bigRes = [];
        p = 1/3;
        n = 18;
        m = 47;
        k = 1000;
        nn = floor(log10(n-1))+1;
        kn = floor(log10(k-1))+1;
        for _ in range(bigK):
                print('-----------------------------------');
                res = [[] for _ in range(n+1)];
                pval_true = [];
                pval_euus = [];
                bin = [nchoosek(n,_)*pow(p,_)*pow(1-p,n-_)*m for _ in range(n+1)];
                print(bin);
                for _ in range(k):
                        sample = list(sum(random()<p for i in range(n)) for i in range(m));
                        frq = freq(sample,n);
                        for i,_ in enumerate(frq):
                                if _:
                                        res[i].append(_);
                        chi2_true = sum((frq[i]-bin[i])**2/bin[i] for i in range(n+1));
                        pval_true.append(gammaincreg(n,chi2_true));
                        chi2_euus = 0;
                        cbin = bin.copy();
                        dir = 0;
                        tn = 0;
                        while frq:
                                cf = 0; cb = 0;
                                while cf < 5 and frq:
                                        cf += frq.pop(dir);
                                        cb += cbin.pop(dir);
                                chi2_euus += (cf-cb)**2/cb;
                                tn += 1;
                                dir = -1-dir;
                        pval_euus.append(gammaincreg(tn-1,chi2_euus));
                ma = max(max(_ or [0]) for _ in res);
                for i,_ in enumerate(res):
                        frq = freq(_,ma);
                        #frq.pop(0);
                        print(('%%%ii'%nn)%i,*[('%%%ii'%kn)%_ if _ else ' '*kn for _ in frq]);
                qval_true = bonferroni(pval_true);
                qval_euus = bonferroni(pval_euus);
                print(qval_euus[:10]);
                print(pval_euus[:10]);
                bigRes.append(sum(_<0.05 for _ in pval_euus));
        print(histf(bigRes));
 

Результат (включая только вывод последнего цикла)

Используется синтаксис Text

-----------------------------------
[0.03180205577614669, 0.2862185019853202, 1.2164286334376107, 3.2438096891669606, 6.08214316718805, 8.515000434063268, 9.224583803568539, 7.906786117344461, 5.4359154556743166, 3.0199530309301754, 1.3589788639185787, 0.49417413233402846, 0.1441341219307583, 0.033261720445559605, 0.005939592936707071, 0.0007919457248942759, 7.424491170883835e-05, 4.367347747578727e-06, 1.2131521521052015e-07]
 0      30
 1     206  40   3   1
 2     362 246  75  24   1   1
 3     103 191 254 191 109  62  28  13   5       1
 4       7  31  84 124 153 187 147 125  74  41  12   9   5
 5           5  13  39  75  88 119 155 153 132  94  62  36  13   8   6   2
 6           5   7  18  45  80 127 136 157 127 111  79  52  30  16   4   5   1
 7           7  20  50  89 124 167 157 130 103  63  48  23   9   6       4
 8      19  61 103 158 189 179 116  79  63  17   7   4   2
 9     146 217 249 161  99  48  25   8   1   1
10     337 258  97  40  10   2   1
11     289  68   7   1
12     122   7   3
13      33
14       5
15
16
17
18
[2.389621366304163e-05, 0.0004849257548949292, 0.002987786017112267, 0.004160743147966028, 0.004160743147966028, 0.02399345014162508, 0.02399345014162508, 0.02420566796471397, 0.02420566796471397, 0.02690354887326138]   
[2.3896213663041628e-08, 9.698515097898584e-07, 8.9633580513368e-06, 1.904709600434731e-05, 2.0803715739830143e-05, 0.0001470063613551422, 0.00016795415099137557, 0.00021379733901034586, 0.00021785101168242572, 0.0002690354887326138]
[[87, 1], [93, 1], [95, 1], [97, 1], [99, 1], [101, 1], [102, 1], [103, 2], [104, 4], [105, 1], [107, 3], [108, 2], [109, 1], [110, 2], [111, 5], [112, 5], [113, 4], [114, 5], [115, 2], [116, 3], [117, 6], [118, 8], [119, 5], [120, 2], [121, 2], [122, 4], [123, 3], [124, 1], [125, 2], [126, 2], [127, 1], [128, 4], [129, 3], [131, 1], [132, 1], [133, 1], [134, 1], [135, 2], [138, 3], [140, 1], [143, 1]]
 

('true' тесты убрал - совсем безобразные результаты дают, всё-таки)

-- 04.07.2024, 23:40 --

EUgeneUS в сообщении #1645205 писал(а):

Сравнить с ожидаемым значением $50$.

Очень странного Вы ожидаете, не находите?

[EUgeneUS]

Geen

                               while cf < 5 and frq:
                                        cf += frq.pop(dir);
                                        cb += cbin.pop(dir);
                                chi2_euus += (cf-cb)**2/cb;

Тут еще ошибка. Цикл должен быть по cb, а не по cf.
Контролируется, чтобы знаменатель не был маленьким.

-- 05.07.2024, 00:57 --

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

[EUgeneUS]

Geen в сообщении #1645206 писал(а):

Очень странного Вы ожидаете, не находите?

Нормального я хочу. P-value так и работает.
1. Строится распределение некой статистики в предположении о верности нулевой гипотезы.
2. Отбрасывается N%-й квантиль.
3. По сути задание критичного уровня P-value - это задание вероятности ошибки первого рода.

Так как обычно эксперимент один, то говорят об уверенности, а не о вероятности. Но у Вас то, как я понял $100 \cdot 1000$ экспериментов.

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

[EUgeneUS]

Geen

EUgeneUS в сообщении #1645203 писал(а):

Завтра, если будет время, сделаю модель на 1000 кейсов в Excel. Посмотрим, насколько будет различаться с Вашим результатом.

Сделал.
Вот результаты типового запуска:
Минимум p-value: $5,69 \cdot 10^{-5}$
Количество, где $\text{p-value} < 0.05$: 46
Количество, где $\text{p-value} < 0.01$: 12

Вот еще:
Минимум p-value: $0.00352$
Количество, где $\text{p-value} < 0.05$: 43
Количество, где $\text{p-value} < 0.01$: 4

И еще
Минимум p-value: $4.66 \cdot 10^{-5}$
Количество, где $\text{p-value} < 0.05$: 53
Количество, где $\text{p-value} < 0.01$: 12

Как видим, отлично работает.
А) никаких $10^{-62}$ нет
Б) количество ошибок первого рода - в полном соответствии с теорией.

Ищите ошибки у себя. На (еще) одну, пока не исправленную, указал выше. Может быть есть ещё.

-- 05.07.2024, 10:39 --

А вот результаты, если не делать группировку

Минимум p-value: $3.14 \cdot 10^{-259}$
Количество, где $\text{p-value} < 0.05$: 92
Количество, где $\text{p-value} < 0.01$: 66

Тут, конечно, так себе работает.

[Geen]

EUgeneUS в сообщении #1645207 писал(а):

Тут еще ошибка. Цикл должен быть по cb, а не по cf.

Да как скажете - это только ухудшит "показатели теста".

EUgeneUS в сообщении #1645227 писал(а):

Ищите ошибки у себя.

Это да, у Вас я их поискать при всём желании не смогу - кроме отдельных цифр и заверений про правильность Вашего пути Вы ничего не предложили....

[EUgeneUS]

Geen в сообщении #1645307 писал(а):

Это да, у Вас я их поискать при всём желании не смогу - кроме отдельных цифр и заверений про правильность Вашего пути Вы ничего не предложили....

1. Вы тут откровенную ерунду написали, мол, я ничего не предложил. Так я и не собирался предлагать.
Я использовал ровно то, что описано в литературе по теме. Ссылку давал.
2. Я Вам по шагам расписал, как, что и с помощью чего делал.
3. Код Вам не предоставил по одной причине - там нет кода. Даже VB.

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

[Geen]

EUgeneUS в сообщении #1645309 писал(а):

Я Вам уже на две ошибки указал

Не две, а одну.

EUgeneUS в сообщении #1645309 писал(а):

Я использовал ровно то, что описано в литературе по теме

Реализацию этого "использования" проверить нет возможности....