Доверительные интервалы полиномиального распределения |
Здравствуйте, гость ( Вход | Регистрация )
Доверительные интервалы полиномиального распределения |
2.04.2020 - 09:24
Сообщение
#1
|
|
Группа: Пользователи Сообщений: 1202 Регистрация: 13.01.2008 Из: Челябинск Пользователь №: 4704 |
Раньше считал (вероятно не совсем корректно) ДИ для долей всегда методами для биномиального распределения. Т.е., например, в ряду абсолютных частот 4-х категорий {13, 35, 18, 7} с суммой n=73 доля первой категории f=13/73*100%=17,8%. Для неё находил 95% ДИ методом Клоппера - Пирсона или в полседнее время методом Джеффриса (байесовский априорный интервал): [10,4; 27,7].
Сейчас решил посчитать ДИ для полиномиального распределения, думал, что раз информации больше, то они Уже будут. Ничего подобного. R-пакет DescTool считает одновременные ДИ для полиномиалного распределения функцией MultinomCI. library(DescTools) x<-c(13,35,18,7) MultinomCI(x) est lwr.ci upr.ci [1,] 0.17808219 0.06849315 0.3006248 [2,] 0.47945205 0.36986301 0.6019947 [3,] 0.24657534 0.13698630 0.3691180 [4,] 0.09589041 0.00000000 0.2184330 По умолчанию считает ДИ методом Сайсона - Глаза по SAS-овскому алгоритму. Всё хуже, чем даже биномиальный Клоппер - Писон, который ругают за консервативность. Видно, что для 7 (9,6%) нижняя граница вообще ноль. Более адекватные результаты даёт только метод Уилсона: > MultinomCI(x, method="wilson") est lwr.ci upr.ci [1,] 0.17808219 0.10713373 0.2812173 [2,] 0.47945205 0.36877454 0.5921840 [3,] 0.24657534 0.16204465 0.3564445 [4,] 0.09589041 0.04722895 0.1849564 Воросы: 1) Каким способом считаете вы? 2) Хочу попробовать сделать бутстреп. Думаю так: многократно пробублировать набор 4 типов в соотношении 13 : 35 : 18 : 7 и извлекать из него с возвратом случайные выборки размером n=73; для каждогго типа потом рассчитать ДИ методом процентилей. Корректно так будет организовать? Сообщение отредактировал nokh - 2.04.2020 - 09:25 |
|
2.04.2020 - 09:49
Сообщение
#2
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Раньше считал (вероятно не совсем корректно) ДИ для долей всегда методами для биномиального распределения. Т.е., например, в ряду абсолютных частот 4-х категорий {13, 35, 18, 7} с суммой n=73 доля первой категории f=13/73*100%=17,8%. Для неё находил 95% ДИ методом Клоппера - Пирсона или в полседнее время методом Джеффриса (байесовский априорный интервал): [10,4; 27,7]. Сейчас решил посчитать ДИ для полиномиального распределения, думал, что раз информации больше, то они Уже будут. Ничего подобного. R-пакет DescTool считает одновременные ДИ для полиномиалного распределения функцией MultinomCI. library(DescTools) x<-c(13,35,18,7) MultinomCI(x) est lwr.ci upr.ci [1,] 0.17808219 0.06849315 0.3006248 [2,] 0.47945205 0.36986301 0.6019947 [3,] 0.24657534 0.13698630 0.3691180 [4,] 0.09589041 0.00000000 0.2184330 По умолчанию считает ДИ методом Сайсона - Глаза по SAS-овскому алгоритму. Всё хуже, чем даже биномиальный Клоппер - Писон, который ругают за консервативность. Видно, что для 7 (9,6%) нижняя граница вообще ноль. Более адекватные результаты даёт только метод Уилсона: > MultinomCI(x, method="wilson") est lwr.ci upr.ci [1,] 0.17808219 0.10713373 0.2812173 [2,] 0.47945205 0.36877454 0.5921840 [3,] 0.24657534 0.16204465 0.3564445 [4,] 0.09589041 0.04722895 0.1849564 Воросы: 1) Каким способом считаете вы? 2) Хочу попробовать сделать бутстреп. Думаю так: многократно пробублировать набор 4 типов в соотношении 13 : 35 : 18 : 7 и извлекать из него с возвратом случайные выборки размером n=73; для каждогго типа потом рассчитать ДИ методом процентилей. Корректно так будет организовать? Для долей надо "восстановить" выборку в виде объектов с признаками замеренными. После этого или 1) перемешиваем признаки (рандомизация) и получаем "доверительный интервал для 0-гипотезы". Или 2) извлекаем "с возвращением" перевыборки и считаем в каждой нужные нам доли, следим за "сходимостью процентилей долей в накопленной выборке", как только достигнутая точность нас устраивает прекращаем. Код > d <- c(rep(1, 13), rep(2, 35), rep(3, 18), rep(4, 7)) > table(d) d 1 2 3 4 13 35 18 7 > length(d) [1] 73 > dd <- replicate(100000, table(sample(factor(d), replace=T))) > str(dd) int [1:4, 1:100000] 12 36 15 10 12 37 16 8 17 40 ... - attr(*, "dimnames")=List of 2 ..$ :8322456 [1:4] "1" "2" "3" "4" ..$ : NULL > sapply(1:nrow(dd), function(i) quantile(dd[i,], probs=c(0.025,0.5,0.975))) [,1] [,2] [,3] [,4] 2.5% 7 27 11 3 50% 13 35 18 7 97.5% 20 43 25 12 PS ну и конечно нарисовать все эти доверительные интервалы можно только в виде тернарных графиков (придется перегруппирование небольшое делать, что бы размерности хватило, хотя можно в виде 3д пирамиды в данном конкретном примере). Сообщение отредактировал p2004r - 2.04.2020 - 10:06 |
|
2.04.2020 - 18:35
Сообщение
#3
|
|
Группа: Пользователи Сообщений: 902 Регистрация: 23.08.2010 Пользователь №: 22694 |
Цитата Каким способом считаете вы? 1. Wilson's CI, которые для данного случая якобы дают более адекватный результат, представляют собой CI для каждого из параметров полиномиального распределения как для биномиального, невзирая на то, что элементарных исходов у эксперимента уже не 2, а 4. Поэтому Уилсоном мы не пользуемся. 2. В оригинальной статье Сайсон* - Глаза (1995) первый из обсуждаемых методов построения доверительного множества для многомерного параметра основан на существовании и оценке такого неотрицательного целочисленного гиперпараметра C, для которого одновременно справедлива система уравнений (6). Поэтому MultinomCI() вовсе не обязаны напоминать о таковых для биномиальных параметров. Резюме: считайте MultinomCI() с установками по умолчанию. 2. Цитата Хочу попробовать ... Корректно так будет организовать? Поскольку мультиномиальное распределение асимметрично, то Эфронов ДИ - это худшее, что можно измыслить в данной ситуации. Патамушта: а) смещение исходной выборки все равно не устраняет; б) непонятно, какие процентили вообще искать: то ли 2.5 - 97.5, то ли 3 - 98, то ли 2 - 97. _____________________________ * Кристина Сайсон - женщина. |
|
4.04.2020 - 05:52
Сообщение
#4
|
|
Группа: Пользователи Сообщений: 1202 Регистрация: 13.01.2008 Из: Челябинск Пользователь №: 4704 |
Благодарю за мнения и код! Попробую всё-таки ещё свой вариант, интересно будет сравнить с результатом р2004r.
По поводу Сайсон - Глаза ничего не читал, но мне решительно не понравился ноль в качестве нижней границы. Получается так: по набору в 73 объекта частота почти 9,5%, а нижняя граница ноль. Причём не 0.0001, что и так нереалистично мало, а вообще 0.00000000. Т.е. по-сути, метод говорит, что несмотря на то, что в выборке у меня оказалось почти 10%, если я продолжу процесс извлечения выборок, то в 95% выборок не обнуружу ни одного объёта такой категории. Не верю. Поэтому более склонен довериться моделированию. Последнее для меня очень затратно по времени написания кодов, но может за самоизоляцию и получится (как ни странно, сейчас времени вообще нет: в НИИ дана команда сидеть дома и писать статьи на год вперёд))), а в универе народ у кого занятий много вообще вешается с этой дистанционкой...) Сообщение отредактировал nokh - 4.04.2020 - 05:54 |
|
4.04.2020 - 05:53
Сообщение
#5
|
|
Группа: Пользователи Сообщений: 1202 Регистрация: 13.01.2008 Из: Челябинск Пользователь №: 4704 |
Не нуда запастил, сорри... Стёр.
Сообщение отредактировал nokh - 4.04.2020 - 05:55 |
|
4.04.2020 - 08:15
Сообщение
#6
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Благодарю за мнения и код! Попробую всё-таки ещё свой вариант, интересно будет сравнить с результатом р2004r. По поводу Сайсон - Глаза ничего не читал, но мне решительно не понравился ноль в качестве нижней границы. Получается так: по набору в 73 объекта частота почти 9,5%, а нижняя граница ноль. Причём не 0.0001, что и так нереалистично мало, а вообще 0.00000000. Т.е. по-сути, метод говорит, что несмотря на то, что в выборке у меня оказалось почти 10%, если я продолжу процесс извлечения выборок, то в 95% выборок не обнуружу ни одного объёта такой категории. Не верю. Поэтому более склонен довериться моделированию. Последнее для меня очень затратно по времени написания кодов, но может за самоизоляцию и получится (как ни странно, сейчас времени вообще нет: в НИИ дана команда сидеть дома и писать статьи на год вперёд))), а в универе народ у кого занятий много вообще вешается с этой дистанционкой...) Возьмите готовые доверительные интервалы отсюда https://cran.r-project.org/web/packages/Ternary/index.html |
|
4.04.2020 - 19:25
Сообщение
#7
|
|
Группа: Пользователи Сообщений: 902 Регистрация: 23.08.2010 Пользователь №: 22694 |
...По поводу Сайсон - Глаза ничего не читал, но мне решительно не понравился ноль в качестве нижней границы. Получается так: по набору в 73 объекта частота почти 9,5%, а нижняя граница ноль. Причём не 0.0001, что и так нереалистично мало, а вообще 0.00000000. Т.е. по-сути, метод говорит, что несмотря на то, что в выборке у меня оказалось почти 10%, если я продолжу процесс извлечения выборок, то в 95% выборок не обнуружу ни одного объёта такой категории. Не верю. Поэтому более склонен довериться моделированию. Все это - следствие Удвоим вектор частот: > x<-2*x > x [1] 26 70 36 14 > MultinomCI(x) est lwr.ci upr.ci [1,] 0.17808219 0.09589041 0.2622103 [2,] 0.47945205 0.39726027 0.5635801 [3,] 0.24657534 0.16438356 0.3307034 [4,] 0.09589041 0.01369863 0.1800185 |
|