Блог

Я пытаюсь найти любой способ для проверки остатков Пуассона, таких как нормали в aov() . В моем гипотетическом примере:

 # For normal distribution
x <- rep(seq(from=10, to=50, by=0.5),6)
y1 <- rnorm(length(x), mean=10, sd=1.5)

#Normality test in aov residuals
y1.av<-aov(y1 ~ x)
shapiro.test(y1.av$res)
#   Shapiro-Wilk normality test
#
#data:  y1.av$res
#W = 0.99782, p-value = 0.7885

  

Звучит глупо, ОК!!

Теперь я хотел бы сделать то же самое приближение, но для распределения Пуассона:

 # For Poisson distribution
x <- rep(seq(from=10, to=50, by=0.5),6)
y2 <- rpois(x, lambda=10)

#Normality test in aov residuals
y2.av<-aov(y2 ~ x)
poisson.test(y2.av$res)
Error in poisson.test(y2.av$res) : 
  'x' must be finite, nonnegative, and integer
  

Есть какой-либо статистический подход для этого?

Спасибо!

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

 #Data
set.seed(123)
# For Poisson distribution
x <- rep(seq(from=10, to=50, by=0.5),6)
y2 <- rpois(x, lambda=10)
  

Теперь мы получаем значения, которые являются элементарными для теста:

 #Values
df <- as.data.frame(table(y2),stringsAsFactors = F)
df$y2 <- as.integer(df$y2)
  

После этого мы должны разделить наблюдаемые значения O и их группы или категории classes . Оба элемента составляют y переменную:

 #Observed values
O <- df$Freq
#Groups
classes <- df$y2
  

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

 #MLE
meanval <- sum(O*classes)/sum(O)
  

Теперь нам нужно получить вероятности каждого класса:

 #Probs
prob <- dpois(classes,meanval)
  

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

 prhs <- 1-sum(prob)
  

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

 #Add probability
prob[length(prob)]<-prob[length(prob)] prhs 
  

С помощью этого мы можем провести тест на соответствие, используя chisq.test() функцию в R . Для этого требуются наблюдаемые значения O и вероятности prob , которые мы вычислили. Просто напоминаю, что этот тест используется для установки неправильных степеней свободы, поэтому мы можем исправить это с помощью формулировки теста, который использует k-q-1 степени. Где k — количество групп и q — количество вычисленных параметров (мы вычислили один параметр с помощью MLE). Далее тест:

 chisq.test(O,p=prob)
  

Результат:

 Chi-squared test for given probabilities

data:  O
X-squared = 7.6692, df = 17, p-value = 0.9731
  

Ключевое значение из теста — это X-squared значение, которое является статистикой теста. Мы можем повторно использовать значение для получения реального p-value (в нашем примере мы имеем k=18 и минус 2, степени свободы равны 16).

p.value Можно получить с помощью следующего кода:

 p.value <- 1-pchisq(7.6692, 16)
  

Результат:

 [1] 0.9581098
  

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