Как оценить пуассоновский процесс с помощью R? (Или: как использовать пакет NHPoisson?)

15

У меня есть база данных событий (то есть переменная дат) и связанных ковариат.

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

Я думаю, что пакет NHPoisson существует только для этой цели; но после 15 часов безуспешных исследований мне все еще далеко не узнать, как им пользоваться.

Черт возьми, я даже пытался читать обе упомянутые книги: Coles, S. (2001). Введение в статистическое моделирование экстремальных значений. Springer. Казелла Г. и Бергер Р.Л. (2002). Статистические выводы. Brooks / Cole.

Один единственный пример в документации к fitPP.fun, кажется, не подходит под мои настройки; У меня нет экстремальных ценностей! У меня просто голые события.

Может ли кто-нибудь помочь мне с простым примером подгонки пуассоновского процесса к параметру с одной ковариатой и предположением, что ? Я заинтересован в оценке и \ alpha . Я предоставляю набор данных из двух столбцов с временами событий (скажем, измеренными в секундах после некоторого произвольного времени t_0 ) и другой столбец со значениями ковариат X ?λИксλзнак равноλ0+αИксλ0αT0Икс

Адам Рычковски
источник
Для тех, кто заинтересован, я работаю над переписыванием этой библиотеки, чтобы улучшить удобство использования. github.com/statwonk/NHPoisson
Statwonk

Ответы:

15

Подгонка стационарного пуассоновского процесса

Прежде всего важно понять, какой тип входных данных нужен NHPoisson.

Прежде всего, NHPoisson нужен список показателей моментов событий. Если мы записываем интервал времени и количество событий в интервале времени, то мы должны перевести его в один столбец дат, возможно, «размывая» даты за интервал, в который они записаны.

Для простоты я предполагаю, что мы используем время, измеренное в секундах, и что «секунда» является естественной единицей .λ

Давайте смоделируем данные для простого стационарного пуассоновского процесса, который имеет событие в минуту:λзнак равно1

lambda=1/60 #1 event per minute
time.span=60*60*24 #24 hours, with time granularity one second

aux<-simNHP.fun(rep(lambda,time.span))

Это simNHP.funделает симуляцию. Мы используем для получения aux$posNHпеременную с индексами моментов срабатывания имитируемого события. Мы можем видеть, что у нас примерно 60 * 24 = 1440 событий, проверяя `length (aux $ posNH).

Теперь давайте перепроектируем с помощью :λfitPP.fun

out<-fitPP.fun(posE=aux$posNH,n=time.span,start=list(b0=0)) # b0=0 is our guess at initial value for optimization, which is internally made with `nlminb` function

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

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

Как только мы это поймем, остальное станет намного проще (по крайней мере, для меня).

Для того, чтобы получить приближение нашего мне нужно взять показатель подогнанной параметра , . И вот приходит еще одна важная часть информации, которая отсутствует в руководстве : с мы вписываемся в логарифм от Х , а не Л себя.λbetaexp(coef(out)[1])NHPoissonλλ

Фиттинг нестационарного пуассоновского процесса

NHPoisson подходит для следующей модели:

λзнак равноехр(пTβ)

пλ

Теперь давайте подготовим нестационарный пуассоновский процесс.

time.span=60*60*24 #24 hours, with time granularity one second
all.seconds<-seq(1,time.span,length.out=time.span)
lambdas=0.05*exp(-0.0001*all.seconds) #we can't model a linear regression with NHPoisson. It must have the form with exp.
aux<-simNHP.fun(lambdas)

Точно так же, как и раньше, aux$posNHмы получим индексы событий, но на этот раз мы заметим, что интенсивность событий экспоненциально уменьшается со временем. И скорость этого уменьшения - это параметр, который мы хотим оценить.

out<-fitPP.fun(tind=TRUE,covariates=cbind(all.seconds),
        posE=aux$posNH,
        start=list(b0=0,b1=0),modSim=TRUE)

Важно отметить, что мы должны ставить all.secondsкак ковариату, а не lambdas. Возведение в степень / логарифмизация выполняется внутри fitPP.fun. Кстати, кроме прогнозируемых значений, функция по умолчанию составляет два графика.

Последняя часть - это функция швейцарского ножа для проверки модели globalval.fun.

aux<-globalval.fun(obFPP=out,lint=2000,
        covariates=cbind(all.seconds),typeI='Disjoint',
        typeRes='Raw',typeResLV='Raw',resqqplot=FALSE)

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

Адам Рычковски
источник
Великолепные объяснения, Адам, большое спасибо. У меня такое впечатление, что мы не можем подобрать модель с двумя группами людей и одной интенсивностью на группу, я прав?
Стефан Лоран
@ StéphaneLaurent Вы можете смоделировать членство в группе как ковариату, и - да, вы можете добавить ковариаты. Могут быть разные интенсивности событий для одной группы и разные для другой. Я никогда не делал ничего подобного.
Адам Рычковски
λя(T)знак равноехр(aя+бT)бя
Адам, может быть, я запутался. Теперь у меня сложилось впечатление, что нет проблем. Я вернусь позже в случае необходимости. Спасибо вам большое за ваше внимание.
Стефан Лоран
поsNЧАС,Nзнак равноTяме,sпaN,беTaзнак равно0)ЕрроряNеяTпп,еUN(поsЕзнак равноaUИксposNH, n = time.span, бета = 0): отсутствует аргумент «start» без значения по умолчанию
vak