Zero-truncated Poisson regression Опции

[ИНО]

Широко используется буржуями для предсказания счетной зависимой переменной, если она в силу дизайна эксперимента не может принимать нулевое значение (ну, например количество прыщей посчитали только на тех больных, которые обратились с жалобами на прыщи). Интерпретация параметров таких моделей заковыриста: https://stats.stackexchange.com/questions/4...ts?noredirect=1

Вопрос: а почему с аналогичной целью нельзя применить обычную регрессию Пуассона, просто предварительно вычтя из всех значений зависимой переменной единицу? Пусть такая модель предсказывает не общее количество счетных единиц, а количество единиц, добавленных к одной обязательно имеющейся. Потом, если надо, добавляем к предсказанной величине единичку, и дело в шляпе! Но, судя по тому, как извращается народ именно с Zero-truncated Poisson regression, и ищет ее программные реализации, очевидно, предложенный мною альтернативный подход неправомерен. Но сколько я не ломал голову, так и не понял, почему. Прошу более головастых подсказать.

И вдогонку еще один маленький вопросик: а как сабж будет грамотно обозвать по-русски? Регрессия Пуассона с усеченным нулем?

[Диагностик]

как сабж будет грамотно обозвать по-русски? Регрессия Пуассона с усеченным нулем?

Ограниченное (усечённое)распределение Пуассона.

Сообщение отредактировал Диагностик - 6.11.2022 - 10:42

[Игорь]

По регрессии Пуассона есть хорошая книга. Свободное скачивание http://qed.econ.queensu.ca/ETM/. С. 475. 2021 года, но по рассматриваемой теме там ничего нет.

Интерпретация параметров таких моделей заковыриста: https://stats.stackexchange.com/questions/4...ts?noredirect=1
... как сабж будет грамотно обозвать по-русски? Регрессия Пуассона с усеченным нулем?

Яндекс.Браузер перевел заглавие ссылки как "Интерпретация коэффициентов регрессии Пуассона с нулевым усечением". По-моему, звучит неплохо.

Сообщение отредактировал Игорь - 6.11.2022 - 17:19

[Игорь]

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

Сообщение отредактировал Игорь - 6.11.2022 - 17:33

[ИНО]

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

А гугль переводить каждый раз по-разному, но все время - как-то криво. "С нулвевым усичением" - не годится, потому как нулевое учечение = отсутствие усечения. А тут фишка в том, что усечение совсем не нулевое, но усекают нули.

Ограниченное (усечённое)распределение Пуассона.

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

Сообщение отредактировал ИНО - 6.11.2022 - 23:54

[Диагностик]

важно отразить, что усечено оно не где-нибудь, а именно в нуле.

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

[ИНО]

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

Сообщение отредактировал ИНО - 7.11.2022 - 01:41

[Диагностик]

Я отвечал на вопрос.

И вдогонку еще один маленький вопросик: а как сабж будет грамотно обозвать по-русски? Регрессия Пуассона с усеченным нулем?

[ИНО]

Ну и? Каким образом данная Вами цитата о том, что распределения бывают усеченные слева, справа и спереди дает ответик на мой конкретный вопросик? Я прекрасно понимаю, что означает "zero-truncated", но не знаю, как это корректно перевести на русский.

[Диагностик]

Усеченное до нуля распределение означает что нуль исключается.

[ИНО]

Спасибо, кэп!

[Диагностик]

Чем богат.

[Игорь]

Хотелось бы высказаться не в ответ, но по теме.
У Johnson с соавт. на 188 описано zero-truncated распределение Пуассона (там же и другие варианты упоминаются). Можно сделать, как у Davidson с соавт. (ссылку см. выше) - рассмотреть схему вычисления zero-truncated регрессии Пуассона: записать ФМП, прологарифмировать, дифференцировать, получить систему нелинейных уравнений, аналитически записать градиент и гессиан, решить систему по Ньютону-Рафсону, оценить качество модели и значимость коэффициентов. Можно пойти и дальше - составить программу на языке Си типа как это реализовано в AtteStat для регрессии Пуассона.
Вопрос, а нужно ли практически. Обосную. Есть зарубежные работы, показывающие принципиально лучшее качество zero-truncated регрессии Пуассона по сравнению с регрессией Пуассона для некоторых наборов данных. Но ситуация похожа на использование ТМФ и критерия Барнарда в анализе таблиц сопряженности - может, и лучше обосновано, но никто не интересуется. Ну еще есть пара примеров из практики, когда эффективные алгоритмы интересуют разве что диссертантов для объема работы. Отсутствие русскоязычных работ (даже официального или хотя бы общепринятого перевода термина не существует, кстати, как и в Барнарде) и скромное количество зарубежных доказывают, что по крайней мере здесь эффективные прикладные алгоритмы анализа данных просто игнорируются научным сообществом, для которого эти алгоритмы по идее должны быть интересны: медики, биологи, экономисты.
Однако так просто оставлять тему нельзя. После Нового года, как обычно, несколько дней нечего делать будет, кроме очистки снега. Сяду на даче в библиотеке, сделаю все сказанное выше хотя бы ради красоты формул и законченности темы. Хоть попрактиковаться. Все-равно преподавание ярких прикладных дисциплин никому здесь не нужно.
Так что спасибо автору темы за идею.

Сообщение отредактировал Игорь - 9.11.2022 - 16:53

[ИНО]

У Johnson с соавт. на 188 описано

Помню, нам преподаватель философии для разрядки обстановки в промежутках между надиктовкой лекции приводил пример того, как давали литературные ссылки в старинных трактатах: "У Платона находим..." - с комментарием от себя: "У какого Платона? Где находим?". По-видимому, Вы приверженец именно той древней школы цитирования

[Игорь]

Как-то в библиотеке заказал редкую статью, типа года 60-го, в которой рассматривалась эквивалентность ранговых критериев. Надеялся, что запросят по МБА. Через какое-то время пригласили придти и выдали листочек с аннотацией. На вопрос, а где же сам текст, ответили, что я оказался первым среди диссертантов, кого не устроила аннотация, а потребовался весь источник. Кстати, о статье - мне потом ее любезно прислал автор, я с ней ознакомился и подарил библиотеке.
Так что - нет, не стиль цитирования, а жизненный опыт. Неужели еще кто-то ходит в библиотеку, сверяется с первоисточником, да просто по ссылкам?

Сообщение отредактировал Игорь - 10.11.2022 - 06:39

[ИНО]

Печальный у Вас жизненный опыт, но не стоит обобщать. По ссылкам ходят и читают, если источник доступен. Если недоступен, ходят, чтобы убедиться в этом. Но первое необходимое условие для этого - наличие этих-самых ссылок в удобоваримом виде. Интересно, сколько есть в природе книг, написанных некоим Ивановым Джонсоном с соавторами с количеством страниц более 188? Прикажете все проверить? И такие с позволения сказать ссылки я встречаю в Ваших постах уже далеко не первый раз.

Сообщение отредактировал ИНО - 10.11.2022 - 07:32

[100$]

Вопрос: а почему с аналогичной целью нельзя применить обычную регрессию Пуассона, просто предварительно вычтя из всех значений зависимой переменной единицу? Пусть такая модель предсказывает не общее количество счетных единиц, а количество единиц, добавленных к одной обязательно имеющейся. Потом, если надо, добавляем к предсказанной величине единичку, и дело в шляпе! Но, судя по тому, как извращается народ именно с Zero-truncated Poisson regression, и ищет ее программные реализации, очевидно, предложенный мною альтернативный подход неправомерен. Но сколько я не ломал голову, так и не понял, почему. Прошу более головастых подсказать.

Надо думать патамушта моделирование суммы X1+X2+,...,+Xn - An независимых одинаково распределенных величин требует перехода к т.н. обобщенному распределению Пуассона (ОПР).
Число n слагаемых в этой сумме является Пуассоновой величиной с параметром lambda.
Как любое безгранично делимое распределение (а ОПР является таковым) каждое обобщенное распределение Пуассона является пределом других ОПР, в т.ч. и "сдвинутых" на некоторую неслучайную алгебраическую величину An. (т.е. характеристическая ф-ция ОПР допускает такой "сдвиг").

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

Я пробовал - для того же набора данных данных получаются существенно другие p, значит, очевидно, есть какой-то теоретический запрет для этого моего ноу-хау.
Понять бы, какой.

с единичной вероятностью.

Сообщение отредактировал 100$ - 12.11.2022 - 18:14

[ИНО]

Благодарю за первый ответ качественный по существу проблемы! Но непонятно: почему мы должны подходить к величине X1 (т. е. первому прыщу из описанного примера) как к случайной, а не как к постоянной? Ведь этот первый прыщ есть примененное условие включения пациента в выборку. Не вижу причин не принимать X1 константой, равной 1, и моделировать только X2.

[100$]

...Но непонятно: почему мы должны подходить к величине X1 (т. е. первому прыщу из описанного примера) как к случайной, а не как к постоянной?

Патамушта в теорвере не определено понятие "постоянная величина".
Ergo ваша единица является случайной. Но, чтобы она постоянно присутствовала в каждом испытании, ее - хочешь, не хочешь - надо объявлять достоверным событием. И всю вероятностную массу распределения сбрасывать в точку 1.

Тогда ваша вероятностная модель получаецца самую малость странной: сначала с единичной вероятностью случается (достоверное) событие - некая "постоянная" величина приняла значение, равное 1, а затем к ней еще с некоторой (пуассоновой) вероятностью выпал 0. Так получается пуассонова 1. Для получения пуассоновой двойки надо 1+0+1, etc. Но, увы вам! Вероятность такого "составного" события равна нулю...

Ведь этот первый прыщ есть примененное условие включения пациента в выборку

Ну, если каждый раз отсутствие признака считать пуассоновым нулем... (типа "зуб даю, этот чувак болен прыщами, хотя на морде у него их ровно 0"), то тогда ваша (неслучайная) выборка условна относительно нуля. Это как раз и есть хорошо нам известное zero-truncated Poisson dist. в чистом виде.

Вообще, вся эта комедия объясняется тем, что описанным вами способом никак нельзя модифицировать носитель распределения (0,1,2,...).
Поэтому ребята вынуждены рассматривать ненужные им нули в рамках условного относительно нуля распределения ("нуль-усеченная модель"), трактовать их как число неудач до первого успеха (отрицательное биномиальное распределение и "negative binomial model"), вводить zero-inflated Poisson regression, вспоминать про изученные Кацем в далеком 1963 году Zero-modified Poisson dist. и нет этому конца.

P.S. А вот, a propos, как бы вы перевели словосочетание Zero-inflated Poisson regression? Настоящий пробный камень начинающего переводчика. Или кошмар.

Сообщение отредактировал 100$ - 13.11.2022 - 16:24

[ИНО]

(типа "зуб даю, этот чувак болен прыщами, хотя на морде у него их ровно 0")

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

Для получения пуассоновой двойки надо 1+0+1, etc. Но, увы вам! Вероятность такого "составного" события равна нулю...

Ту что-то совсем не понял. Двойка у нас никакая не пуассоновская, она получается сложением постоянной 1 с пуассоновской 1. Предположением, что мы не заворачиваем никого в дверях, но за первый "прыщик" на лице пациента считаем... его нос. Он есть у всех и каждого, в отличие от второго и последующего прыщиков, которые действительно прыщики. Врач - инопланетянин, причем очень тупой, и не умеет отличать нос от прыщиков, а может только посчитать общее количество выпуклостей на лице. Вполне себе попадутся ему и с одной, и с двумя и с тремя и т. д., никаких нулевых вероятностей, за исключением собственно нулевых знаний, которых гарантировано не будет (ну, разве то герой Гоголя на прием придет ). По-моему вполне себе адекватная вероятностная модель выходит.

[100$]

Ту что-то совсем не понял. Двойка у нас никакая не пуассоновская, она получается сложением постоянной 1 с пуассоновской 1.

Ну вот, началось в колхозе утро: я ему про вероятности, он мне - про школьную алгебру. Так, чего доброго, и до таблицы умножения дойдет.
А какая же она? Детерминированная, чтоле?

Двойка - самая что ни на есть пуассонова:
- во-первых, несомненно случайная, потому как зависит от случайного второго слагаемого;
- во-вторых, много ли вы знаете дискретных распределений, носителем которых является {0,1,2,3...}?

Если она не пуассоновская, то какая? Из распределения Скеллама, что ли?
А, может быть, Накагами? Трейси-Видома? Уишарта? Уточните, когда не лень.

P.S. Вы распределение-то сгенерируйте: result<-rpois(1000,.685). При параметре распределения Lambda<1 0 - самое вероятное значение. Как от него избавляться? Распределение Пуассона - это же все-таки закон редких событий.

[ИНО]

Что-то я Вас с каждым постом все меньше понимаю. Зачем избавляться от пуассоновского от 0? Он ничуть не мешает, даже наоборот необходим! 1(константа) + 0 (пуассоновский) =1 - как раз самое часто значение в моих подсчетах "прыщиков". И частоты всех последующих значений (1+1=2, 1+2=3...) вполне согласуются с гипотезой о добавлении постоянной единицы к случайной пуассоновской величине. Этот первый прыщик - непременное условие для осмотра больного, достоверное событие, а вот уже количество дополнительных прыщиков - случайная величина. Наезд на алгебру тоже не понял. Прибавление постоянного свободного члена в регрессионной модели Вас тоже смущает?

Сообщение отредактировал ИНО - 14.11.2022 - 10:07

[100$]

Что-то я Вас с каждым постом все меньше понимаю. Зачем избавляться от пуассоновского от 0? Он ничуть не мешает, даже наоборот необходим! 1(константа) + 0 (пуассоновский) =1 - как раз самое часто значение в моих подсчетах "прыщиков". И частоты всех последующих значений (1+1=2, 1+2=3...) вполне согласуются с гипотезой о добавлении постоянной единицы к случайной пуассоновской величине. Этот первый прыщик - непременное условие для осмотра больного, достоверное событие, а вот уже количество дополнительных прыщиков - случайная величина.

Ну, а что тут непонятного? Вы сконструировали некоторое дискретное распределение. Для дискретного распределения можно явно указать вероятность любого события. Вот и укажите вероятность появления суммы вида 1(константа) + 0 (пуассоновский) =1.
Пока этого не сделано, я буду считать это школьной математикой.

Наезд на алгебру тоже не понял. Прибавление постоянного свободного члена в регрессионной модели Вас тоже смущает?

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

Сообщение отредактировал 100$ - 14.11.2022 - 13:04

[ИНО]

Что такое "школьная математика" я не в курсе, равно как и о причине Вашего уничижительного отношения к алгебре.

Вот и укажите вероятность появления суммы вида 1(константа) + 0 (пуассоновский) =1

Та же, что и вероятность 0 в обычной пуассоновской регрессионной модели, построенной по той же выборке, из каждой варианты которой вычтена 1. Формулой это не напишу, увы.

Нет. Потому что в регрессионной модели константа - это такой же (оцениваемый) регрессионный коэффициент, и так же точно является случайной величиной

Может быть оценен, а может быть и взят с неба из теоретических предпосылок, и коэффициенты, кстати, тоже, так что оцениваться будет только ошибка. И это все равно будет статистическая модель.

И да тут все ж форум о прикладной статистике, а не о математической (это разные дисциплины, хотя и связанные). О последней пусть д. ф.-м. н. семиэтажными формулами дискутируют, теоретически обосновывая очередной элегантный чудо-метод имени себя (то, что на практике он может толком не работать, потому что "забыли про овраги" - это уже дело десятое ) . Так что объяснение "на пальцах" для простых смертных типа меня всяко приветствуется.

Сообщение отредактировал ИНО - 15.11.2022 - 09:09

[100$]

Та же, что и вероятность 0 в обычной пуассоновской регрессионной модели, построенной по той же выборке, из каждой варианты которой вычтена 1. Формулой это не напишу, увы.

Ну, если на пальцах, то вы сконструировали некое условное относительно неслучайной величины (единицы) распределение и утверждаете, что в нем условная вероятность равна безусловной.
Хотя это ниоткуда не следует.
Вот, например, если мы захотим сдвинуть стандартное нормальное распределение, мы же не будем городить такую словесную конструкцию "а вот мы дождемся с вероятностью 1 единицу, а потом посмотрим, как реализуется ст.норм.расп." Мы этот сдвиг просто добавим к параметру распределения. И оно будет симметрично теперь уже относительно 1.
Вот бы и с Пуассоном также точно.

[ИНО]

Как обычно, я попробовал развеять теоретический сумрак моделированием. Решил: какая из двух моделей будет иметь более точный прогноз на смоделированных данных, той буду пользоваться на практике. Все ж я биолог, а не математик. Но возник конфуз: прогнозы у сабжа и "моего безграмотного ноу-хау" и их ошибки оказались практически неотличимыми. Более того, при лямбда<1 (а это больше похоже на реальные данные) моя модель работает даже чуточку лучше! Хотя возможно тут дело не непосредственно в лямбде, а в том, что мой алгоритм моделирования усекает объемы выборок, сгенерированных при малых лямбда, но тогда это означает, что расово правильная модель хуже моей ведет себя на малых выборках, так что хрен редьки не слаще.

Уменьшение n на описанную картину не влияет (правда, при n меньше 100 и очень малых лямбда, алгоритм сталкивается с невозможностью оценки модели из-за полного отсутствия ненулевых значений и стопорится, так что в этом диапазоне условий исследование провести не удалось. Наверное, по-хорошему надо применять вместо цикла сопли-функции или как их там, а я в этом не силен. Но для понимания проблемы и имеющееся сойдет.

Итак, прошу указать на мою ошибку, не в теории, которую как тут выяснилось, я постичь не в силах, а именно в моделировании!

Код

# генерация синтетических данных:
lambda<-seq(0.01, 2, 0.01)
simmat<-matrix(NA, 200, 1000)
for(i in 1:200){
lambda_i<-lambda[i]
simmat[i,]<-rpois(1000, lambda_i)}

# предложенная мной раскритикованная выше модель:
prediction<-rep(NA, 200)
MAE<-rep(NA, 200)
for(i in 1:200){
x<-simmat[i,]
df<-as.data.frame(x)
df<-subset(df, x>0)
df$modx<-df$x-1
mod<-glm(modx~1, data=df, family="poisson")
pred<-predict(mod, type="response")+1
prediction[i]<-pred
MAE[i]<-mean(abs(pred-df$x))}

# модель c расово правильно усеченным нулем:
library(VGAM)
prediction2<-rep(NA, 200)
MAE2<-rep(NA, 200)
for(i in 1:200){
x<-simmat[i,]
df<-as.data.frame(x)
df<-subset(df, x>0)
mod<-vglm(x~1, data=df, family="pospoisson")
pred<-predict(mod, type="response")
prediction2[i]<-pred
MAE2[i]<-mean(abs(pred-df$x))}

# пргнозы обеих моделей (моя - красненькая):
plot(lambda, prediction, type="n")
lines(lambda, prediction, col="red")
lines(lambda, prediction2, col="blue")

# средние абсолютные ошибки обеих моделей:
plot(lambda, MAE, type="n")
lines(lambda, MAE, col="red")
lines(lambda, MAE2, col="blue")

# резульаты их сравнения при различных лямбда:
plot(lambda, MAE>MAE2)

Сообщение отредактировал ИНО - 18.11.2022 - 10:49

Эскизы прикрепленных изображений

 

[Диагностик]

при лямбда<1

А что такое лямбда?

[ИНО]

Параметр распределения Пуассона. В моем случае имеется в виду распределение, из которого были сгенерированы синтетические данные. Следует отметить, что варианты с нулевым значением исключены из последующего анализа. По идее, именно исходное распределение, их включающее, пытается реконструировать zero-truncated Poisson regression. Мой же подход работает по-другому: в нем оценивается обычное распределение Пуассона, по выборочным значениям, из которых была вычтена единица. У этого распределения уже иная лямбда, которая здесь не приводится. Да, для простоты в обоих случаях построена регрессия к 1, т, е. просто произведена оценка параметров распределения. Регрессионный анализ как таковой - это уже следующий шаг.

[Диагностик]

Параметр распределения Пуассона.

То есть матожидание. А как здесь ТЕХ использовать?

Сообщение отредактировал Диагностик - 18.11.2022 - 14:32

[100$]

Наконец-то наш естествоиспытатель сам признался, что переливать из пустого в порожнее ему сподручнее безо всякой теории )

Однако теория гласит, что ежели моделировать Пуассоном счетные данные, то с необходимостью надо принять, что эти данные измерены в абсолютной шкале. А в ней, как это знает не только взрослый, но даже карапуз, допустимы только тождественные преобразования, т.е. преобразования вида f(x)=x, переводящие каждый элемент шкалы в самого себя. Ни тебе сложения, ни умножения, ни возведения в степень, ни извлечение квадратного корня. Как из таких данных можно вычесть 1 непонятно.

Далее. В Пуассоново распределение мутирует биномиальное распределение при малой вероятности успеха и неограниченном возрастании N. А что есть биномиальное распределение?
Вот выпало в результате эксперимента из 100 бросков монеты 49 орлов и 51 решка. Как из этих данных вычесть единицу? Это ж фальсификацыя ).

А модель... Ну, что модель: у нее же нет "защиты от дурака". Если она не столкнулась с вычислительными проблемами, то она чего-то там посчитает. Ну и что?

[comisora]

2 ИНО

У меня не хватает компетенции ответить на Ваш исходный вопрос (есть более квалифицированные форумчане). Однако посчитал возможным дать ссылки на материалы, которые мне попались при поиске информации по вопросу:

1. Difference between shifted distribution and zero-truncated distribution - по ссылке демонстрация разницы;
2. Simulation of the Shifted Poisson Distribution with an Application to the CEV Model - теоретическая статья, в которой обсуждается генерация такого распределения.
3. Will adding a constant to a random variable change its distribution? - обсуждение того, как меняется распределение при добавлении константы.

[comisora]

2 ИНО

Ещё немного материала.

"In this chapter, we discuss models for zero-truncated and zero-inflated count data. Zero truncated means the response variable cannot have a value of 0. A typical example from the medical literature is the duration patients are in hospital. For ecological data, think of response variables like the time a whale is at the surface before re-submerging, counts of fin rays on fish (e.g. used for stock identification), dolphin group size, age of an animal in years or months, or the number of days that carcasses of road-killed animals (amphibians, owls, birds, snakes, carnivores, small mammals, etc.) remain on the road. These are all examples for which the response variable cannot take a value of 0.

On their own, zero-truncated data are not necessarily a problem. It is the underlying assumption of Poisson and negative binomial distributions that may cause a problem as these distributions allow zeros within their range of possible values. If the mean is small, and the response variable does not contain zeros, then the estimated parameters and standard errors obtained by GLM may be biased. In Section 11.2, we introduce zero-truncated Poisson and zero-truncated negative binomial models as a solution for this problem. If the mean of the response variable is relatively large, ignoring the truncation problem, then applying a Poisson or negative binomial (NB) generalised linear model (GLM), is unlikely to cause a problem. In such cases, the estimated parameters and standard errors obtained by Poisson GLM and truncated Poisson GLM tend to be similar (the same holds for the negative binomial models).

In ecological research, you need to search very hard to find zero-truncated data. Most count data are zero inflated. This means that the response variable contains more zeros than expected, based on the Poisson or negative binomial distribution. A simple histogram or frequency plot with a large spike at zero gives and early warning of possible zero inflation. This is illustrated by the graph in Fig. 11.1, which shows the numbers of parasites for the cod dataset that was used in Chapter 10 to illustrate logistic regression. In addition to presence and absence of parasites in cod, Hemmingsen et al. (2005) also counted the number of parasites, expressed as intensity" (Zuur AF, Ieno EN, Walker N, Saveliev AA, Smith GM. Mixed effects models and extensions in ecology with R [Internet]. New York, NY: Springer New York; 2009 [cited 2022 May 30]. (Statistics for Biology and Health). Available from: http://link.springer.com/10.1007/978-0-387-87458-6).

Дальнейший поиск следует вести по следующим ключевым словам: "the Shifted Poisson distribution", "the Positive Poisson distribution". Так по второму ключевому слову нашёл краткое упоминание в работе автора библиотеки VGAM (Yee TW. Vector Generalized Linear and Additive Models [Internet]. New York, NY: Springer New York; 2015 [cited 2020 Oct 13]. (Springer Series in Statistics). Available from: http://link.springer.com/10.1007/978-1-4939-2818-7).

Далее приведу пример, демонстрирующий разницу между распределениями:

CODE

set.seed(32167)
table(VGAMdata::rpospois(1000,0.5))

1 2 3 4
763 194 39 4

set.seed(32167)
table(rpois(1000, 0.5))

0 1 2 3
601 303 78 18

Рекомендую к просмотру презентацию, в которой объяснены некоторые идеи.

[ИНО]

100$, я уже видел Ваши 100500 теоретических аргументов в пользу того, что шмель летать моя модель работать не может, спасибо. 100501-й и 100502-й уже погоды не сделают. Теперь я отошел от изначальной формулировки вопроса и интересуюсь иным: почему результаты моего моделирования все же показывает, что шмель летает моя модель работает не хуже (а местами даже самую малость лучше лучше), чем аэроплан zero-truncated Пуассон. Прошу заметить, что смоделированные мною данные являются классическим случаем zero-truncated Пуассона, который обязан адекватно обрабатываться соотвествующей функцией из пакета VGAM (если, конечно, ее не дурак писал, но, вроде, не жалуются) и поэтому Ваши отсылки к условиям применимости идут лесом. Имеется факт: на синтетическом примере моя "неправильная" модель справляется не хуже общепризнанной для смоделированной ситуации расово правильной. Но это лишь в том случае, если я не ошибся в моделировании. в чем совсем не уверен. Поэтому прошу проверить мой код и указать на ошибку, если таковая имеется. Теоретико-философских же аргументов довольно, я их услышал, а далее уже становится непонятно, кто из нас на самом деле тут из пустого в порожнее переливает. Вот я, например, код написал и даже анимацию сотряпал, а Вы?

[ИНО]

comisоra, в очередной раз спасибо за полезные ссылки! Судя по ним, то что я здесь предложил, - это именно shifted Poisson distribution. И я не первый, кто задавался вопросом о возможности его применения для "положительной регрессии": https://stats.stackexchange.com/questions/2...hifting-vs-trun. Но, к сожалению, в этой теме обсуждение сместилось в другую область и заданные автором вопросы так и остались без ответа. Увы, больше ничего полезного по запросу "shifted Poisson distribution" найти не смог. Отличия zero-truncated от обычного Poisson distribution я, в принципе представляю, и это не то, о мне хотелось бы сейчас почитать.

Итак, решил я еще немного "из пустого в порожнее попереливать", а именно немного модифицировать код из своего позапрошлого поста, так чтобы на каждом шаге модели строились по одинаковому числу "наблюдений". В данном примере оно установлено на 10, что приближает условия учений к боевым . За одно пришлось увеличить объем изначальной пуассоновской выборки до 2000, так как из 1000 при малых лямбда было, что не набиралось и десятка единиц. Таким образом мы полностью исключили влияние эффекта неравновеликости выборок на поведение моделей и можем утверждать, что моя неправильная модель таки обходит "правильную" именно при малых лямбда, да и при больших не отстает. В целом моя модель почти в три раза чаще имеет меньшую среднюю абсолютную погрешность, чем "правильная".

Код

# генерация синтетических данных:
lambda<-seq(0.01, 2, 0.01)
simmat<-matrix(NA, 200, 10)
for(i in 1:200){
lambda_i<-lambda[i]
sim<-rpois(2000, lambda_i)
sim<-sim[sim>0]
simmat[i,]<-sample(sim, size=10)
}

# предложенная мной раскритикованная выше модель:
prediction<-rep(NA, 200)
MAE<-rep(NA, 200)
for(i in 1:200){
x<-simmat[i,]
df<-as.data.frame(x)
df$modx<-df$x-1
mod<-glm(modx~1, data=df, family="poisson")
pred<-predict(mod, type="response")+1
prediction[i]<-pred
MAE[i]<-mean(abs(pred-df$x))}

# модель c расово правильно усеченным нулем:
library(VGAM)
prediction2<-rep(NA, 200)
MAE2<-rep(NA, 200)
for(i in 1:200){
x<-simmat[i,]
df<-as.data.frame(x)
mod<-vglm(x~1, data=df, family="pospoisson")
pred<-predict(mod, type="response")
prediction2[i]<-pred
MAE2[i]<-mean(abs(pred-df$x))}

y<-MAE>MAE2
length(y[y==TRUE]) # число случаев, в которых расово правильная модель оказалась лучше моей
length(y[y==FALSE] ) # число случаи, в которых моя модель оказалась лучше расово правильной

Сообщение отредактировал ИНО - 19.11.2022 - 16:41

[Диагностик]

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

С чем сравнивались обе модели?

[ИНО]

Друг с другом.

[Диагностик]

И какая из них считается правильной?

[ИНО]

Этот вопрос лучше адресовать $100, поскольку мне , в отличие от него, неведомы критерии истинно правильности моделей. Могу лишь сказать, что по результатам тестирования на синтетических данных у моей модели оказалась чуть ниже ошибка предсказания. Но это очень малое отличие, на практике совершенно несущественное. Непонятно иное: почему вопреки теоретическим соображениям, моя модель не оказалась значительно хуже, чем та, которая, по идее, должна идеально описывать именно то распределение, данные из которого я генерирвал

Сообщение отредактировал ИНО - 20.11.2022 - 12:08

[comisora]

Непонятно иное: почему вопреки теоретическим соображениям, моя модель не оказалась значительно хуже, чем та, которая, по идее, должна идеально описывать именно то распределение, данные из которого я генерирвал

Код

sim<-rpois(2000, lambda_i)

2 ИНО

Может попробуете погенерировать VGAMdata::rpospois()?

Мне представляется некорректным Ваш подход исходного вычитания единицы в одном случае, так как в другом случае Вы этого не делаете. Считаю, что проверять модели надо в одинаковых условиях, а именно после генерации (два варианта):
1) Без фильтрования нулей добавить ко всему ряду +1;
2) После фильтрации нулей работать с тем, что получилось без добавлений и вычитаний.
Также попробуйте сравнить модели по другим характеристикам (AIC, MSE и т.п.).

Сообщение отредактировал comisora - 20.11.2022 - 12:49

[100$]

Я на том стою, что базовое распределение Пуассона, из коего командой

Код

sim<-sim[sim>0]

тупо повыбрасывали все нули, не имеет никакого отношения к Shift'у и не имитирует вероятностных свойств случайной величины,
"которая по условиям эксперимента не может принимать значение 0" (с).
Кроме того, такая команда равносильна заявлению, что автор пакета VGAM - дурак, не понимающий простых вещей: он-то зачем-то заморочился генерацией ZTP-распределения (rpospois()), а тут всего-то и делов: из обычного распределения Пуассона знай себе нули выбрасывай - и дело с концом.

Вам comisora вроде бы ясно указал на разницу этих распределений, ан не в коня корм.

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

Изначально речь шла о том, что для ZTP распределения мы сначала прогоняем ZPT-модель, получаем "заковыристо интерпретируемые" коэффициенты, а потом на этом же распределении пытаемся "обмануть" обычную регрессию Пуассона.

Сообщение отредактировал 100$ - 20.11.2022 - 14:03

[ИНО]

100$, а как по-вашему получается этот "позитивный Пуассон" в реальности? По-моему, именно отбрасыванием всех нулевых значений (они просто не попадают в нашу, как однажды выразился 100$, неслучайную выборку) - ровно то же, что делает мой генератор. Конечно, на самом деле реальное изначальное распределение может и не быть пуассоновским (скорее всего, в большинстве реальных случаев это будет "zero-inflated") , но для простоты пока что условимся считать его таковым. В чем согласен: к "шифту" мой генератор действительно не имеет отношения, но чудесным образом на нем шифт-модель работает лучше, чем ZTP.

comisora, вечерком попробую этот rpospois, хотя покамест слабо представляю, чего он такого может делать, принципиально отличного от моего ручного усекновения нулей, и почему это что-то должно лучше описывать реальный процесс порождения данных. И, хотя у меня есть некоторые сомнения относительно преимущества квадратичных ошибок при оценке счетной величины, на всякий случай посчитаю также и MSE, и RMSE. А вот AIC тут уж точно лишний - некого штрафовать.

[100$]

100$, а как по-вашему получается этот "позитивный Пуассон" в реальности? По-моему, именно отбрасыванием всех нулевых значений (они просто не попадают в нашу, как однажды выразился 100$, неслучайную выборку) - ровно то же, что делает мой генератор. Конечно, на самом деле реальное изначальное распределение может и не быть пуассоновским (скорее всего, в большинстве реальных случаев это будет "zero-inflated") , но для простоты пока что условимся считать его таковым. В чем согласен: к "шифту" мой генератор действительно не имеет отношения, но чудесным образом на нем шифт-модель работает лучше, чем ZTP.

Вы как-то очень по-своему понимаете разницу между "с.в. не может принимать нулевые значения" и "может, но мы не будем считать это исходом эксперимента и все равно потом повыкидаем".

Поэтому галиматья получается уже на этапе генерации этой вашей "синтетики": я прошу вас сгенерировать мне выборку объемом 2000 из ZTP. Вы честно генерируете выборку из распределения Пуассона (нуачотакова?) объемом 2000, в которой оказывается порядка 1200 нулей. Вы их выбрасываете, и с честными глазами предъявляете мне выборку объемом 800. А я просил 2000...
Словом сделаете так раз-другой... "А потом ваши рыжие кудри примелькаются, и вас просто начнут бить".

P.S. Ваша клоунада несколько затянулась.

Сообщение отредактировал 100$ - 20.11.2022 - 15:21

[ИНО]

Слышал, что люди применительно к уровню знаний о конкретном предмете делятся на категории:

1) знает, что знает;
2) думает, что знает;
3) думает, что не знает;
4) знает, что не знает;
5) не знает, что знает;
6) не знает, что не знает

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

100$, с каждой вашей новой нападкой на меня Вы становитесь все менее похожим на умного человека. Поясняю: умный человек на вашем месте, прежде чем искрометно блеснуть своей гениальностью и посрамить оппонента, на всякий пожарный проверил бы, а действительно ли по результатам работы мой генератор "распределения-уродца", существенно отличается от расово правильного генератора ZTP из VGАM (автор которого, как выяснилось, точно не дурак), и, если да, то действительно ли эти различия влияют на эффективность обсуждаемых моделей. А то вдруг окажется, что генератор вообще не причем? Всего-то пара примитивных строчек кода и меньше минуты вычислений - а какая замечательная страховка от конфуза! Но Ваше ЧСВ, очевидно, настолько велико, что всякий раз дает железобетонную априорную уверенность в собственной правоте, не нуждающуюся в каких-либо плебейских проверках. Оттого раз за разом Вы попадаете в, скажем так, неловкие ситуации (хотя, возможно, с высоты Вашего персонального олимпа они и не кажутся таковыми, но окружающие, уж поверьте, каждый промах Акелы замечают и мотают на ус). Предыдущая была с некой функцией CWtest() Вашего авторства. Поэтому, если претендуете на высокое звание человека, принадлежащего к первой из вышеозначенных категорий, придется такие досадные промахи исключить полностью путем дотошных перепроверок любого собственного тезиса, позиционируемого как непреложная истина, до его выноса на публику. Моя-то позиция в этом плане несравненно проще - это категория ?4 Так что могу себе позволить в высказываниях заблуждаться безо всяких ограничений, желая быть поправленным более знающими людьми. Увы, покамест в данной теме это мое желание удовлетворено не было.

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

Только диаграммку приложу, что такое x и y - сами догадайтесь, я в Вас верю!

Сообщение отредактировал ИНО - 20.11.2022 - 18:01

Эскизы прикрепленных изображений

 

[ИНО]

1) Без фильтрования нулей добавить ко всему ряду +1;
2) После фильтрации нулей работать с тем, что получилось без добавлений и вычитаний.

Можно много чего посчитать, но нужно понимание целей процесса. Цель предложенных Вами методов неясна, какие реальные ситуации Вы предлогаете при помощи них смоделировать? Первый предложенный Вами подход как раз и есть генерация shifted Poisson distribution, на статью о которой вы дали ссылку в предыдущем посте. Кстати, я так и не понял, зачем его автор предлагал аж три разных метода генерации и почему не было достаточно просто прибавить или отнять константу, т. е. собственно произвести этот самый shift. Но на практике мне не встречались явления, для которых такое распределение можно было бы допустить в качестве процесса порождения данных. Реальные данные, с которыми я работаю как раз и есть результат "фильтрации нулей", таков естественный процесс их сбора. Именно его я попытался восседать в своем генераторе, приняв допущение о том, что до отбраковки распределение было пуассоновским. И именно это фактически делает ZTP (пусть и немножко иным путем).

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

Проверил и снова удивился: обычная пуассоновская модель на данных, сгенерированных rpospois(), вопреки теоретическим соображениям показала себя тоже неплохо, только на больших выборках (n=1000) она систематически отстает от остальных двух (моей сдвинутой и расово правильной ZTP), на малых же (n=20) положение трех моделей на пьедестале почета меняется случным образом в зависимости от конкретной выдачи датчика случайных чисел, но моя неправильная лидирует чаще других. В любом случае разница мизерная.

Читал я, что разница между прогнозами обычной модели Пуассона и модели ZTP тем больше, чем меньше лямбда. Мол при больших лямбда можно смело применять обычную пуассоновскую модель, забив на отсутствие нулей, а вот при малых - ни-ни. По теоретическим соображениям логично. Но проверка на синтетических данных показала странное. На приведенных ниже диаграммах попарно сравниваются средние абсолютные ошибки трех моделей: моей сдвинутой (MAE), ZTP (MAE2) и обычной пуассоновской (MAE3). Данные сгенерированы расово правильной функцией rpospois(), n=1000. На удивление обычная пуассоновская модель обошла остальные именно на малых лямбда, где по идее должна быть неприменима от слова "совсем". Но обратите внимание на порядок величин разностей, этот мизер может быть в большей степени связан с программной реализацией, чем реальными свойствами распределений. Графики прогнозов и ошибок практически идентичны, так что не привожу.

Сообщение отредактировал ИНО - 20.11.2022 - 20:01

Эскизы прикрепленных изображений

 

[100$]

А что не так с CWTest()'ом?
Там, вроде бы, перестановочный тест Фишера-Питмана не мог отвергнуть гипотезу менее чем на 5%-ном уровне.

Все остальные тесты отвергали бы ее за счет чумовой ошибки I рода.

P.S. Кстати, а что это за Большой Брат, который тут за мной наблюдает и чего-то там себе на что-то наматывает? Вы там всей палатой сгрудились, штоле?

[comisora]

Цель предложенных Вами методов неясна...

Чтобы зависимые переменные были одинаковые. В таком случае сравниваться будет только результат влияния параметра family.

Реальные данные, с которыми я работаю как раз и есть результат "фильтрации нулей", таков естественный процесс их сбора. Именно его я попытался восседать в своем генераторе, приняв допущение о том, что до отбраковки распределение было пуассоновским.

Раз в естественном процессе сбора есть нули в результате отбраковки, почему бы не использовать hurdle-модель?

[ИНО]

100$, А Вы тему просмотрите (благо тут из свежих негусто) и уведите, что не так с CWTest()'ом. "Чумовая ошибка" первого рода была только у него, у всех прочих - второго, что и следовало ожидать при столь малых размерах выборок. И только это не имеющее аналогов в мире изобретение отличлось! Причем, пока я не продемонстрировал данный конфуз моделированием, Вы уж чуть было не написали пакет имени себя, даже не удосужившись проверить работоспособность собственного кода (а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков). Теперь, похоже, наблюдается рецидив. Вместо бестолковых рассуждений о Большом брате и палатах, лучше б R запустили. и убедились что мой обруганный Вами генератор ровно то же, что и хваленый rpospois() - подсказка относительно рисунка из поста #43 для тех, кто сам не догадался.

comisora, эта гибрид из логичстической модели и ZPT, применяется. наоборот, когда в выборке нулей овер дофига. Наверное, можно расковырять готовую функцию, и оторвать половину гибридной модели, но получим все тот же pospoisson.

Зависимые переменные в моем коде и так одинаковые, столь огорчившие теоретика арифметические фокусы с единицей - это уже часть регрессионной модели. Кстати, собственно к регрессии мы покамест еще перешли, этот шаг запланирован на завтра. Может, с появлением независимых переменных модели наконец-то поведут себя по-разному, и неправильные облажаются. Хотя в теории загвоздка должна была заключаться именно в той части процесса, которая уже проделана, - в моделировании распределения остатков - однако покамест теория что-то плохо контачит с практикой, так что я уже не удивлюсь любым сюрпризам.

Сообщение отредактировал ИНО - 20.11.2022 - 22:44

[comisora]

comisora, эта гибрид из логичстической модели и ZPT, применяется. наоборот, когда в выборке нулей овер дофига. Наверное, можно расковырять готовую функцию, и оторвать половину гибридной модели, но получим все тот же pospoisson.

Я не вижу причин отказываться от hurdle-модели даже когда мало нулей, но изначально их быть не должно. Такая составная модель, как мне представляется, может справиться с процессом отбраковки (1 или 0), а далее Пуассоном с усечённым нулём моделировать оставшиеся данные.

CODE

n = 10
lambda = 0.5

set.seed(32167)
x <- qpois(p = runif(n, min = dpois(0, lambda)), lambda) #rpospois()

set.seed(32167)
y <- qpois(p = runif(n, min = 0), lambda)

set.seed(32167)
z <- rpois(n, lambda)

print(x);print(y);print(z)

[x] 1 1 1 1 1 2 2 1 1 2
[y] 0 0 0 0 0 1 2 0 0 2
[z] 0 0 0 0 0 1 2 0 0 2

Я пока не придумал как адаптировать Ваш код для доказательства/опровержения утверждения, что Ваш генератор и rpospois() даёт одно и тоже. Пока только демонстрирую разницу между rpospois() и rpois() с невозможностью арифметическими операциями свести одно распределение к другому. Даже если мы возьмём и создадим последовательность Пуассона без нулей, то она не всё равно не совпадёт при большом количестве наблюдений и lambda менее 10. Дальнейшее увеличение lambda приводит к расхождению рядов генерации.

CODE

n = 10000
lambda = 9
set.seed(32167)
x <- qpois(p = runif(n, min = dpois(0, lambda)), lambda)
set.seed(32167)
y <- qpois(p = runif(n, min = 0), lambda)
set.seed(32167)
z <- rpois(n, lambda)
table(x-y);table(x-z);table(y-z)

[100$]

100$, А Вы тему просмотрите (благо тут из свежих негусто) и уведите, что не так с CWTest()'ом. "Чумовая ошибка" первого рода была только у него, у всех прочих - второго, что и следовало ожидать при столь малых размерах выборок. И только это не имеющее аналогов в мире изобретение отличлось! Причем, пока я не продемонстрировал данный конфуз моделированием, Вы уж чуть было не написали пакет имени себя, даже не удосужившись проверить работоспособность собственного кода (а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков). Теперь, похоже, наблюдается рецидив. Вместо бестолковых рассуждений о Большом брате и палатах, лучше б R запустили. и убедились что мой обруганный Вами генератор ровно то же, что и хваленый rpospois() - подсказка относительно рисунка из поста #43 для тех, кто сам не догадался.

Ба! Затянувшаяся клоунада сменилась затянувшейся истерикой.

Так, все-таки, кто это там мои "промахи" себе на что-то наматывает? От чьего имени вы тут вещаете? На гельминтоз давно проверялись?

Вам, видимо, невдомек, что гарантированно ошибку I рода удерживают только перестановочные критерии, вычисляющие точное распределение. У всех остальных (с аппроксимациями) она будет либо чумовая, либо чуть менее чем. Но это общее место в статистике. Зачем проверять теорему Пифагора линейкой?
Так что никакого "конфуза".

Что касается "пакета имени себя", то мотивация была немного иная.
Вот ссылка, прочитав которую я усомнился в своей способности проделать то, о чем там написано (н-р, п.п. 6-7). Если бы тема получила продолжение, то объяснил бы, в чем дело. Ну, вот объяснил теперь. Лучше поздно, чем никогда.

И, кстати, что такого нерабочего в моем коде?
Что вы разумеете под "работоспособностью"?
Что это у вас там так подгорело?
Здесь же не форум программистов: здесь бесконечное число юзеров способны написать что угодно на чем угодно быстрее и лучше меня. И чо?

Зачем оно нужно вообще и мне в частности? Интересно сначала решить задачу, а потом кропать программы.

...И только это не имеющее аналогов в мире изобретение отличлось!
...(а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков)...

А вот идиотничать я вам не позволю.


Сообщение отредактировал 100$ - 21.11.2022 - 00:45

[ИНО]

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

Что касается KWTeest(), приведу Вашу же цитату:

Я тут чуток помонтекарлил для равных выборочных объемов (n1=n2) от 3 до 100 способность CWTest() удерживать номинальный уровень значимости (alpha=,05).
Получается, что сколько-нибудь осмысленное его применение возможно только при n1=n2>60, п.ч. сходимость к предельному распределению очень медленная. Ребята, это же ужасно.

Пусть каждый решает для себя, нужен ли ему такой ущербный критерий при наличие давно известных намного лучших альтернатив. Как по мне, это просто мусор. Но не надо втирать, будто иные критерии ведут себя схожим образом, ибо это наглое вранье, CWTest() - воистину аналоговнет. Ладно, предположим, что великий 100$ не стал создавать пакет имени себя не потому, что осознал ущербность этого своего критерия, а чисто в силу неспособности освоить техническую сторону процесса (вернее "сомнения в собственных силах", ведь, как известно, истинный мэтр может все, но имеет право иногда подвергать эту всемогучесть легкому сомнению ). Но почему же он в этой теме не написал ни строчки кода R в подтверждение своих обвинительных заявлений? Неужто и в этой способности через полгода возникли сомнения? Или просто лень одолела?

Сообщение отредактировал ИНО - 21.11.2022 - 10:44

[ИНО]

Я не вижу причин отказываться от hurdle-модели даже когда мало нулей, но изначально их быть не должно. Такая составная модель, как мне представляется, может справиться с процессом отбраковки (1 или 0), а далее Пуассоном с усечённым нулём моделировать оставшиеся данные.

Пожалуйста, проиллюстрируйте примером. А то я сосем не понимаю, что вы имеете в виду.

Я пока не придумал как адаптировать Ваш код для доказательства/опровержения утверждения, что Ваш генератор и rpospois() даёт одно и тоже.

Это просто:

# генерация синтетических данных:
lambda<-seq(0.01, 2, 0.01)
x<-rep(NA, 200)
y<-rep(NA, 200)
for(i in 1:200){
lambda_i<-lambda[i]
x[i]<-mean(rpospois(1000, lambda_i))
sim<-rpois(1000000, lambda_i)
sim<-sim[sim>0]
y[i]<-mean(sample(sim, size=1000))
}

plot(x, y)
abline(0, 1, col="red")

Здесь x - выдача rpospois(), а y - моего генератора. Фактически он делает ровно то же, но способом, который в более явном виде имитирует реальный процесс порождения "позитивнопуассоновских" данных. Конечно, полного соответствия не будет, т. к. у меня датчик псевдослучайных чисел работает в двух разных функциях и я не знаю как в такой ситуации корректно задать одинаковое семя обоим генераторам. Но и полученной диаграммы вполне достаточно для осознания их эквивалентности. Поскольку мы имеем дело с однопараметрическим распределением, достаточно оценки одного параметра, но для пущей убедительности можно проделать то же самое с любой другой, скажем заменить mean() на sd() - картина не изменится.

Пока только демонстрирую разницу между rpospois() и rpois() с невозможностью арифметическими операциями свести одно распределение к другому. Даже если мы возьмём и создадим последовательность Пуассона без нулей, то она не всё равно не совпадёт при большом количестве наблюдений и lambda менее 10. Дальнейшее увеличение lambda приводит к расхождению рядов генерации.

Это изначально понятно. Именно поэтому я и удивляюсь: почему оценивая через сдвинутое пуассоновское распределение данные, сгенерированные усеченным пуассоновским распределением (т. е. как бы некорректно) я получаю практические те же оценки (а то и чуть более точные), как если бы действовал по феншую, выбрав заведомо верную спецификацию модели. Это действительно странно.

[Диагностик]

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

Для какого лямбда?

[ИНО]

Для любого в интервале от 0,001 до 2 (на больших не проверял, ибо неинтересно).

[ИНО]

Итак, переходим к собственно регрессии. Пока что пусть будет один бинарный предиктор. В этот раз для генерации данных взяты лямбды, похожие на оцененные в одной из моих реальных задач, и объем выборки того же порядка. В этот раз мы будем сравнивать модели на 1000 разных реализаций распределений с одними и теми же параметрами. Для разнообразия вместо MAE взята MSE, но и с MAE картина ровно та же - модели практически не различаются по точности прогнозов. Традиционно под номером 1 идет моя оригинальная шифт-модель, под номером 2 - ZTP, под номером 3 - обычная пуассоновская

Код

pr1<-rep(1, 50)
pr2<-rep(2, 50)
pr<-as.factor(c(pr1, pr2))

TestPoisReg<-function(pr)
{
resp1<-rpospois(50, 0.6)
resp2<-rpospois(50, 1.35)
resp<-c(resp1, resp2)
df<-data.frame(pr, resp)
df$modresp<-df$resp-1
mod1<-glm(modresp~pr, data=df, family="poisson")
prediction1<-predict(mod1, type="response")+1
MSE1<-mean((prediction1-df$resp)^2)
mod2<-vglm(resp~pr, data=df, family="pospoisson")
prediction2<-predict(mod2, type="response")
MSE2<-mean((prediction2-df$resp)^2)
mod3<-glm(resp~pr, data=df, family="poisson")
prediction3<-predict(mod3, type="response")
MSE3<-mean((prediction1-df$resp)^2)
return(c(MSE1, MSE2, MSE3))
}
x<-replicate(1000,  TestPoisReg(pr=pr))
MSE1<-x[1,]
MSE2<-x[2,]
MSE3<-x[3,]

boxplot(MSE1, MSE2, MSE3)

Сообщение отредактировал ИНО - 21.11.2022 - 16:38

Эскизы прикрепленных изображений

 

[ИНО]

А вот в p разница существенная:

Код

pr1<-rep(1, 50)
pr2<-rep(2, 50)
pr<-as.factor(c(pr1, pr2))

TestPoisReg2<-function(pr)
{
resp1<-rpospois(50, 0.6)
resp2<-rpospois(50, 1.35)
resp<-c(resp1, resp2)
df<-data.frame(pr, resp)
df$modresp<-df$resp-1
mod1<-glm(modresp~pr, data=df, family="poisson")
prediction1<-predict(mod1, type="response")+1
p1<-coef(summary(mod1))[2,4]
mod2<-vglm(resp~pr, data=df, family="pospoisson")
p2<-coef(summary(mod2))[2,4]
mod3<-glm(resp~pr, data=df, family="poisson")
p3<-coef(summary(mod3))[2,4]
return(c(p1, p2, p3))
}
x<-replicate(1000,  TestPoisReg2(pr=pr))
p1<-x[1,]
p2<-x[2,]
p3<-x[3,]

boxplot(p1, p2, p3)

Ну, третью модель (обычный Пуассон) можно отправлять сразу в мусорное ведро, потому что она пытается предсказать нули, которых нет, построенные по ней ДИ будут невозможными, соответственно и p тоже должно быть неправильное. А вот как сделать выбор между ?1 и ?2, понять какая из них больше врет в данном случае, я покамест не сообразил. Буду благодарен за подсказку. И все ж странно, почему при таких серьезных различиях в объясняющей составляющей (которая была ожидаема), в составляющей прогностической все три модели оказались практически одинаковыми.

Сообщение отредактировал ИНО - 21.11.2022 - 18:33

Эскизы прикрепленных изображений

 

[ИНО]

За отсутствием альтернативных предложений решил попробовать ответить на последний поставленный мною вопрос самостоятельно, а именно при помощи моего оригинального метода, предложенного здесь: http://forum.disser.ru/index.php?showtopic...ost&p=27378

Теорема имени меня : чем ближе распределение p-значений какого-либо статистического критерия при справедливости нулевой гипотезы к равномерному, тем ближе номинальный уровень значbмости к реальному. При этом смещение моды распределения влево от значения 0,5 сигнализирует о завышенной ошибке первого рода, а вправо - о завышенной ошибке второго рода. Доказывает (или опровергает) пусть $100 - у него голова большая он шибко умный

Исследуем при помощи предложенного метода три вышеупомянутые модели:

Код

pr1<-rep(1, 50)
pr2<-rep(2, 50)
pr<-as.factor(c(pr1, pr2))

TestPoisReg2<-function(pr)
{
resp1<-rpospois(50, 1)
resp2<-rpospois(50, 1)
resp<-c(resp1, resp2)
df<-data.frame(pr, resp)
df$modresp<-df$resp-1
mod1<-glm(modresp~pr, data=df, family="poisson")
prediction1<-predict(mod1, type="response")+1
p1<-coef(summary(mod1))[2,4]
mod2<-vglm(resp~pr, data=df, family="pospoisson")
p2<-coef(summary(mod2))[2,4]
mod3<-glm(resp~pr, data=df, family="poisson")
p3<-coef(summary(mod3))[2,4]
return(c(p1, p2, p3))
}
x<-replicate(1000,  TestPoisReg2(pr=pr))
p1<-x[1,]
p2<-x[2,]
p3<-x[3,]

hist(p1)
abline(v=0.5, col="red")

hist(p2)
abline(v=0.5, col="red")

hist(p3)
abline(v=0.5, col="red")

Ожидаемо у модели ?3 (обычный Пуассон) с достигаемым уровнем значимости для эффекта предиктора все плохо - он сильно завышен. В этот раз, когда предиктор является мусорным, оно даже хорошо, однако и для по-настоящему влиятельного предиктора он тоже будет завышен, что уже совсем нехорошо. Так что использовать обычную регрессию Пуассона для данных с отсутствующими нулями с целью оценки влияния независимых переменных однозначно недопустимо (хотя с задачей прогноза в такой ситуации она вполне справляется, что весьма странно, но все результаты моделирования покамест говорят именно об этом). А вот с моделями ?1 (моя шифтонутая) и ?2 (ZTP) все не так однозначно: в зависимости от n и лямбда, задаваемых генератору, распределение p красивее то у одной, то у другой из них, но разница всегда небольшая. В случае, представленном на прикрепленных гистограммах лидирует расово правильная усеченная модель, а моя дает малость заниженное p, но в некоторых иных, мною смоделированных, напротив, моя модель имеет почти идеальное распределение, а расово правильная демонстрирует признаки небольшого завышения. Может поэкспериментировать самостоятельно.

Сообщение отредактировал ИНО - 23.11.2022 - 17:52

Эскизы прикрепленных изображений

 

[100$]

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

[ИНО]

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

Почему Вы упомянуты? Так Вы тут единственный (оставшийся), кто полагает (ИМХО необоснованно), что в совершенстве понимает теорию математической статистики. Так что можете считать это обращение знаком особого уважения к вашему интеллекту. Но нет, так нет.

И, раз уж начали к никам придираться (распоследний аргумент в любой интернет-дискуссии), на свой посмотрите. Это что, тариф за час вашего умственного труда? Лучше быть трехбуквенным, чем ценником





Сообщение отредактировал ИНО - 23.11.2022 - 20:31

[Игорь]

Врачи-аспиранты, для которых создан форум, ознакомившись с показанными выше результатами исследований, могут составить впечатление о непостижимости темы для медика или биолога. Сугубо личное мнение, просто такие исследования не нужно обсуждать здесь. А нужно [в рассматриваемом случае] представлять совет в виде:
1. Есть данные, являющиеся результатами практических измерений в эксперименте.
2. Смотрим, какая статистическая модель (модели) подходит (подходят) для их описания (с точностью до коэффициентов).
3. Идентифицируем модель - вычисляем оптимальные значения коэффициентов, при необходимости значимость параметров, включенных в модель.
4. Проверяем соответствие модели экспериментальным данным. Часто применяется критерий AIC. Если подходит несколько моделей, выбираем лучшую.
5. Используем построенную модель для прогнозирования и в публикациях. Всё красиво, все довольны.

Есть еще одна идея. Ранее, когда компании-разработчики статистического программного обеспечения (ПО) рассылали демонстрационные материалы (диски, брошюры) по почте, получил от StaXact очень интересные картонные карточки формата меньше A4 (наверное, B5) с голубой каемкой. На лицевой стороне карточки была описана проблема - реальный случай из медицинской практики, включая ссылки. На оборотной стороне - результаты статистического исследования с интерпретацией результатов. Идея понравилась - найди пример и сделай так же.

Ведь смысл научного медико-билогического исследования не в развитии математико-статистических методов, а во всестороннем анализе объектов исследования, развитии методик исследований и корректном применении в качестве инструментов готовых математических алгоритмов и ПО, как и приборов, оборудования. Если хотим разрабатывать математические алгоритмы, кодировать ПО и создавать приборы - для обсуждения этого есть другие соответствующие площадки.

Не нужно пугать аспирантов.

[Диагностик]

Для любого в интервале от 0,001 до 2 (на больших не проверял, ибо неинтересно).

Распределения близки, можно выбирать любое. На графиках теоретические вероятности.



Сообщение отредактировал Диагностик - 24.11.2022 - 10:42

[ИНО]

Диагностик, при других лямбдах разница будет заметнее. См. первую ссылку из этого поста comisora: http://forum.disser.ru/index.php?s=&sh...ost&p=27854

Здесь не обсуждается отличия в форме между сдвинутым и усеченным распределениями Пуассона, оно самоочевидно. Обсуждается разница в работоспособности (настаиваю на этом термине) регрессионных моделей, основанных на том и другом. Результаты моделирования (если я не напутал чего, но покамест никто не поправил) свидетельствуют, что в задачах оценивания разница может быть существенной, а вот в задачах прогноза ее практически нет (более того, можно использовать даже обычную модель Пуассона). Последнее странно, и я все ж надеюсь, что напортачил в коде, ибо иначе придется таки уверовать в "непостижимость темы для биолога".

Игорь, ИМХО Ваш совет "не исследовать" звучит странно. Эдак, слепо следуя рекомендациям всяких там гуру (не будем называть ники) можно таких дел наворотить, что ой. Все ж считаю, что любой исследователь должен понимать хотя бы на начальном уровне суть каждого метода сбора и обработки данных, который он применяет. Обратное чревато неверными выводами, что в сфере медицинских исследований особенно опасно. И если аспирант - настолько нежная и чувствительная личность, что прочтение этой темы может его напугать, то занят он явно не своим делом. Такому лучше попробовать реализовать себя в искусствах, а не в науке (да, я в курсе, что медицина - отчасти тоже искусство, но это уж точно искусство не для пугливых, речь об иных). Из тонких экзальтированных особ получаются отличные поэты, художники и композиторы, но никак не ученые.

Что же касается биологов, в особенности экологов, то их роль в развитии современной прикладной статистики трудно переоценить. Созданные ими методы (нередко для решения одной конкретной специфической задачи, с которой не справились ранее существующие) зачастую в последствие оказываются намного более полезными и универсальными в практическом применении, чем всякие там равномерно наиболее мощные сферокони, заботливо выращиваемые теоретиками от статистики.

[Диагностик]

Обсуждается разница в работоспособности (настаиваю на этом термине) регрессионных моделей, основанных на том и другом. Результаты моделирования (если я не напутал чего, но покамест никто не поправил) свидетельствуют, что в задачах оценивания разница может быть существенной,

Обе будут работать. Нужно проверять на фактическом статматериале.

[Диагностик]

Игорь, ИМХО Ваш совет "не исследовать" звучит странно.

Я бы вас поддержал, если бы ваши личные неуклюжие попытки исследования не выглядели так не профессионально.

В случае, представленном на прикрепленных гистограммах лидирует расово правильная усеченная модель

Предлогаю остановиться на этом.


Сообщение отредактировал Диагностик - 25.11.2022 - 06:54

[ИНО]

Остановиться в развитии - весьма плохое предложение.

Диагностик, что позволено Юпитеру, не позволено быку. Если 100$, несмотря излишнюю обидчивость и воинственность в отстаивании собственных ошибочных суждений (как было здесь наглядно доказано для тех, кто понимает), все еще заслуживает того, чтобы быть услышанным, поскольку его знания о теории математической статистики достаточно высоки (значительно выше моих, хотя и ниже того непогрешимого уровня, который он сам сам назначил), то Вы - обычный дилетант с поверхностными знаниями, от которого я не припомню хотя бы одного глубокомысленного и полезного поста (хотя я зарегистрирован на данном форуме и недавно, но почитываю его уж давненько). Поэтому уж кому-кому, а не Вам судить об уровне моих исследований. Если мои выводы ошибочны, что я вполне допускаю (и едва ли не в этом посте об этом напоминаю!), то укажите на конкретные ошибки (причем не высокими философствованиями о свойствах случайных величин, понять которые способны не только лишь все), а конкретными строчками кода. Моя проблема проста: данные моделирования противоречат теории. Принимая в качестве постулата справедливость теории, я прошу элементарного: укажите мои ошибки в моделировании. 100$ указал (на мой самописный генератор усеченного распределения Пуассона), но попал пальцем в небо, поскольку проверка показала, что этот генератор не хуже расово правильного вгамовского (хоть и медленнее). Уверен, было б у него иное объяснение, уже б давно гарцевал тут на лихом коне, подняв его на знамя. Но вместо этого тихонько отсиживается в сторонке, лишь изредка высыпая порцию оскорблений в адрес моей личности. Это есть не что иное как примитивнейшая из уловок ведения дискуссии, применение которой оппонентом однозначно сигнализирует о том, что аргументы по существу у него закончились. Так что если он не сумеет напрячь свою гениальность и таки родить состоятельное объяснение обсуждаемого феномена, остается лишь констатировать его фиаско. А куда уж Вам с вашими привычными сентенциями капитана Очевидность? Оттого ваши обвинения пусты и смехотворны. Я всегда стараюсь не переходить на личности в научной дискуссии, но, ежели оппонент сделал это первым, то сам виноват - что посеял, то и пожнет.

Вообще, к глубочайшему сожалению, приходится констатировать, что единственный оставшийся активный пользователь данной ветки, который действительно стремится от чистого сердца давать полезные советы, при этом проводя громадную работу с литературой, это comisora. За что ему нижайший поклон. Прочие, увы, либо перестали посещать форум, либо забронзовели. Ддаже те, кто не имел ни малейшего основания для этого, типа Вас, Диагностик.

Сообщение отредактировал ИНО - 25.11.2022 - 11:39

[Диагностик]

объяснение обсуждаемого феномена

Какого феномена?

[100$]

Я бы вас поддержал, если бы ваши личные неуклюжие попытки исследования не выглядели так не профессионально.

Я худею, дорогая редакция: на арене появился второй коверный и сцепился с первым.
Диагностик, только 6.38, а вы уже нетрезвы?
Впрочем, есть более лестная для вас версия: если учесть что в 2 часа ночи вы (возможно, мимоходом) еще присутствовали на форуме, то тогда все более или менее понятно: чего не напишешь спросонок...

P.S. А ведь я предупреждал, что два юзера типа "Диагностик" для этого форума многовато. Так сказать "Боливару не снести двоих"...

[Диагностик]

Профессор на семинаре замечает, что один из студентов не совсем трезв. Говорит ему:
- Так, студент, вы пьяны! Немедленно покиньте семинар!
- Да что вы, профессор?! Я совершенно трезв!
- Ах, трезв?! Тогда напишите мне формулу вычисления лапласиана в тороидальных координатах.
Студент выходит к доске и пишет длинную-длинную формулу. Профессор долго сверяется со справочником и убеждается, что всё написано правильно.
- Ну что ж, я впечатлён! Я эту формулу и на трезвую-то голову не помню.
Студент: "Так на трезвую и я не помню!"

Сообщение отредактировал Диагностик - 25.11.2022 - 14:07

[100$]

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

[Диагностик]

P.S. А ведь я предупреждал, что два юзера типа "Диагностик" для этого форума многовато.

Но не прислушались. Вот теперь и полюбуйтесь к чему это привело, прозорливый вы наш.

А ведь я предупреждал

Сообщение отредактировал Диагностик - 25.11.2022 - 14:36

[100$]

Но не прислушались. Вот теперь и полюбуйтесь к чему это привело, прозорливый вы наш.

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

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

[Диагностик]

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

После первого его сообщения я также его поблагодарил в ЛС и сообщил, что его информация должна снять вашу проблему.

Сообщение отредактировал Диагностик - 25.11.2022 - 14:54

[100$]

Уверен, было б у него иное объяснение, уже б давно гарцевал тут на лихом коне, подняв его на знамя. Но вместо этого тихонько отсиживается в сторонке, лишь изредка высыпая порцию оскорблений в адрес моей личности. Это есть не что иное как примитивнейшая из уловок ведения дискуссии, применение которой оппонентом однозначно сигнализирует о том, что аргументы по существу у него закончились. Так что если он не сумеет напрячь свою гениальность и таки родить состоятельное объяснение обсуждаемого феномена, остается лишь констатировать его фиаско.

Товарищ! Добрее надо быть: вам рано еще здесь кого-либо понукать.
Так что не нукай: не запрягал.

Если, конечно, вы не рассчитываете на то, что я по щелчку ваших пальцев удалюсь от зла, и со с страшной силой начну доить козла...

[Диагностик]

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

А что вы понимаете под

...задачах прогноза..

[Диагностик]

Моя проблема проста: данные моделирования противоречат теории.

Вы можете это предметно показать? Какие гипотезы и на каком уровне вами проверялись?

[ИНО]

Специально два дня не читал эту тему (хотя видел, что новые ответы в ней есть) в наивной надежде, что за это время напишут хоть что-то умное. Увы, надежда не оправдалась: написанное в мое отсутствие оказалось еще более жестяной жестью, чем до того. За сим не вижу никакой пользы в дальнейшем вникании в ваши разборки в стиле "Горца" за честь остаться последним оставшимся единственно правильным гуру. Тем более, что они весьма унылы: ни тебе фехтования на мечах, ни аргументированной научной дискуссии, ни строчки кода, а одно только пустословное хамство. Если б я жаждал беседы подобного "научного уровня", сходил бы вечерком к ближайшей наливайке: несмотря на все ваши старания, тамошних завсегдатаев в этом специфическом, с позволения сказать, искусстве аргументации вам ни в жизнь не переплюнуть. Жив буду - вернусь. Когда-нибудь. Постарайтесь к тому моменту не сжить друг друга со свету с форума. Не хотелось бы, чтобы к моему возвращению на ветке остался только один Макклауд - тогда она станет совсем уж печальной и бесполезной.

Сообщение отредактировал ИНО - 27.11.2022 - 16:39

[100$]

Специально два дня не читал эту тему (хотя видел, что новые ответы в ней есть) в наивной надежде, что за это время напишут хоть что-то умное. Увы, надежда не оправдалась: написанное в мое отсутствие оказалось еще более жестяной жестью, чем до того. За сим не вижу никакой пользы в дальнейшем вникании в ваши разборки в стиле "Горца" за честь остаться последним оставшимся единственно правильным гуру. Тем более, что они весьма унылы: ни тебе фехтования на мечах, ни аргументированной научной дискуссии, ни строчки кода, а одно только пустословное хамство. Если б я жаждал беседы подобного "научного уровня", сходил бы вечерком к ближайшей наливайке: несмотря на все ваши старания, тамошних завсегдатаев в этом специфическом, с позволения сказать, искусстве аргументации вам ни в жизнь не переплюнуть. Жив буду - вернусь. Когда-нибудь. Постарайтесь к тому моменту не сжить друг друга со свету с форума. Не хотелось бы, чтобы к моему возвращению на ветке остался только один Макклауд - тогда она станет совсем уж печальной и бесполезной.

Я всего лишь (вежливо) указал Диагностику на недопустимость подобных речеваний в ваш адрес.
Диагностик (еще вежливее) указал на допустимость подобных словесных эскапад. Можно сказать, оставил за собой право.

Какие еще "фехтования на мечах"? Какие "научные дискуссии"? Какие "строчки кода"? Причем здесь старые галоши?
Ничего не понимаю.

Но ясно одно: ваша логика - того-с: демонстративно вам изменяет.
Вам надо поскорее сменить рюмочную.
В этой "наливайке" вам постоянно наливают что-то не то.



Сообщение отредактировал 100$ - 29.11.2022 - 14:24

[Диагностик]

Топикстартёру мною лично был задан вопрос. Он его либо не понял, либо не услышал, либо проигнорировал ("Не сказать ещё хужей.")

[100$]

Тут вот еще какая штука выяснилась: на неотрицательном носителе k=1,2,3,... определено также логарифмическое распределение. Ссылка
Интересно, чтобы его сгенерировать, тоже надо у Пуассоновского распределения тупо нули отсечь?
Вот вопрос вопросов.

[Диагностик]

См. Генераторы дискретно распределенных случайных величин, Логарифмическое распределение. https://habr.com/ru/post/265321/

[Игорь]

См. Генераторы дискретно распределенных случайных величин, Логарифмическое распределение. https://habr.com/ru/post/265321/

Спасибо. Красивая библиотека на правильном языке программирования. Кодом, конечно, пользоваться не буду, но документацию посмотрю.

Сообщение отредактировал Игорь - 1.12.2022 - 06:50

[Игорь]

Говорил выше, что библиотека примеров для исследователей - неплохая идея ...

Свежий корректный пример регрессии Пуассона https://nplus1.ru/news/2022/12/02/from-bike-to-transplant

[100$]

Я продолжу эту тему поиском ответа на основной вопрос ТС: "можно ли вычесть единицу из zero-truncated Poisson distribution и прогнать обычную регрессию Пуассона?".

Ссылка на то, что это вообще такое -ZTP вот

Под нумером [6] в перечне ссылок фигурирует страница с кодами/способами генерации ZTP distribution. Ссылка

Там рассмотрены два способа генерации - (A) и (B)

Способ (A) "Rejection sampling":

Код

n<-1000; T<-3.5
  Y<-rpois(n,T); Y0<-Y[Y>0]; r<-(n - length(Y0))

  while(r>0)
{ Y<-rpois(r,T); Y0<-c(Y0,Y[Y>0]); r<-(n - length(Y0)) }

ясно указывает на то, что "генерация выборки с отклонением" представляет собой не однократное усечение нулей, как доказывал ТС, а "итерационный" возобновляемый процесс многократного усечения, причем параметр r на каждой новой генерации/итерации зависит от количества нулей, отброшенных на предыдущей итерации.

Так я ответил на свой собственный вопрос о том, как же обеспечить заявленный объем выборки в 2000, если из нее выбросить ~1200 нулей.
А то ТС с ответом затруднился.

Способ генерации (B)

n<-1000 # desired size of sample
T<-3.5 # pre-truncation mean of Poisson
U<-runif(n) # the uniform sample
t = -log(1 - U*(1 - exp(-T))) # the "first" event-times
T1<-(T - t) # the set of (T-t)

X <- rpois(n,T1)+1 # the final truncated Poisson sample

явно указывает на то, что при вычитании единицы из ZTP получится выборка, каждый элемент которой сгенерирован с параметром T1=(T-t), что означает разную распределенность элементов такой выборки. Ergo считать ее выборкой независимых одинаково распределенных пуассоновских величин нет оснований. И она не описывается единым параметром распределения, который в регрессии выступает в качестве зависимой (моделируемой) переменной.
По-видимому, это и есть тот самый теоретический "запрет", превращающий ZTP в распределение - "мутант". Просто бессмыслица получится. Хотя на выборках объемом 10 это может и не бросаться в глаза.

Сообщение отредактировал 100$ - 5.12.2022 - 00:03