Как я могу оценить плотность нулевого параметра в R?

10

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

set.seed(1)
x <- c(rlnorm(100),rep(0,50))
hist(x,probability=TRUE,breaks = 25)

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

lines(density(x), col = 'grey')

Есть density(... from, to)аргументы, но они, похоже, только усекают вычисления, а не изменяют окно, так что плотность в 0 согласуется с данными, как видно на следующем графике:

lines(density(x, from = 0), col = 'black')

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

Есть ли альтернативы этой функции, которые обеспечили бы лучший расчет плотности в нуле?

введите описание изображения здесь

Abe
источник

Ответы:

14

Плотность бесконечна в нуле, потому что это включает дискретный шип. Вам необходимо оценить пик, используя пропорцию нулей, а затем оценить положительную часть плотности, предполагая, что она гладкая. KDE вызовет проблемы на левой стороне, потому что это придаст вес отрицательным значениям. Одним из полезных подходов является преобразование в журналы, оценка плотности с использованием KDE, а затем преобразование обратно. См. Wand, Marron & Ruppert (JASA 1991) для справки.

Следующая функция R выполнит преобразованную плотность:

logdensity <- function (x, bw = "SJ") 
{
    y <- log(x)
    g <- density(y, bw = bw, n = 1001)
    xgrid <- exp(g$x)
    g$y <- c(0, g$y/xgrid)
    g$x <- c(0, xgrid)
    return(g)
}

Тогда следующий даст сюжет, который вы хотите:

set.seed(1)
x <- c(rlnorm(100),rep(0,50))
hist(x,probability=TRUE,breaks = 25)
fit <- logdensity(x[x>0]) # Only take density of positive part
lines(fit$x,fit$y*mean(x>0),col="red") # Scale density by proportion positive
abline(v=0,col="blue") # Add spike at zero.

введите описание изображения здесь

Роб Хиндман
источник
Спасибо за ваш ответ, но я в замешательстве - вы говорите: «Оцените пик, используя пропорцию нулей», но постройте его без границ. имеет ли шип дискретную высоту или он бесконечен, если он дискретен, это ? P(X=0)
Абэ
Это смесь дискретного распределения и непрерывного распределения. Если график представлен плотностью, то шип бесконечен (на самом деле дельта-функция Дирака). Иногда люди изображают дискретную часть как функцию вероятностной массы (тогда шип имеет высоту ), а непрерывную часть - как функцию плотности. Это, вероятно, делает изображение лучше, но оно включает в себя две разные шкалы. P(X=0)
Роб Хиндман
это пригодится. к вашему сведению: похоже, что, хотя bw = "SJ" влияет на плотность в нетрансформированном пространстве, плотность записи одинакова при использовании "SJ" и значения по умолчанию "nrd0" ... Я собираюсь прочитать ссылку на SJ: "Sheather and Джонс (1991). Надежный метод выбора полосы пропускания на основе данных для оценки плотности ядра ". jstor.org/stable/2345597
Abe
4

Я согласен с Робом Хиндманом, что вам нужно разбираться с нулями отдельно. Существует несколько методов оценки плотности ядра переменной с ограниченной поддержкой, в том числе «отражение», «перенормировка» и «линейная комбинация». Похоже, что они не были реализованы в densityфункции R , но доступны в пакете Бенна Янна kdensдля Stata .

универсальный
источник
1

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

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

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

Грег Сноу
источник
0

Вы можете попробовать уменьшить пропускную способность (синяя линия для adjust=0.5), введите описание изображения здесь

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


источник
Есть ли другой метод, который вы бы порекомендовали?
Абэ
@Abe Ну, это зависит от того, что вы хотите сделать ...