Подгонка кривой плотности к гистограмме в R

91

Есть ли в R функция, которая соответствует кривой гистограмме?

Допустим, у вас была следующая гистограмма

hist(c(rep(65, times=5), rep(25, times=5), rep(35, times=10), rep(45, times=4)))

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

Это довольно простой вопрос, но я не могу найти ответ для R в Интернете.

r histogram curve-fitting r-faq user5243421
источник
Вы хотите найти такие m и s, чтобы гауссово распределение N (m, s) соответствовало вашим данным?
SteinNorheim,
Я не уверен, что это значит ...> _>
user5243421
10
@mathee: Я думаю, он имел в виду m = среднее и s = стандартное отклонение. Гауссово распределение - это еще одно название нормального распределения.
Питер Мортенсен,

Ответы:

154

Если я правильно понял ваш вопрос, то вам, вероятно, понадобится оценка плотности вместе с гистограммой:

X <- c(rep(65, times=5), rep(25, times=5), rep(35, times=10), rep(45, times=4))
hist(X, prob=TRUE)            # prob=TRUE for probabilities not counts
lines(density(X))             # add a density estimate with defaults
lines(density(X, adjust=2), lty="dotted")   # add another "smoother" density

Отредактируйте спустя некоторое время:

Вот немного более нарядная версия:

X <- c(rep(65, times=5), rep(25, times=5), rep(35, times=10), rep(45, times=4))
hist(X, prob=TRUE, col="grey")# prob=TRUE for probabilities not counts
lines(density(X), col="blue", lwd=2) # add a density estimate with defaults
lines(density(X, adjust=2), lty="dotted", col="darkgreen", lwd=2)

вместе с графиком, который он производит:

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

Дирк Эддельбюттель
источник
3
+1 - можно ли сделать наоборот, т.е. настроить график плотности в соответствии с гистограммой?
vonjd
2
Я предлагаю указать дополнительный параметр, lines(density(X,na.rm= TRUE)поскольку вектор может содержать значения NA.
Anirudh
30

Это просто с ggplot2

library(ggplot2)
dataset <- data.frame(X = c(rep(65, times=5), rep(25, times=5), 
                            rep(35, times=10), rep(45, times=4)))
ggplot(dataset, aes(x = X)) + 
  geom_histogram(aes(y = ..density..)) + 
  geom_density()

или имитировать результат решения Дирка

ggplot(dataset, aes(x = X)) + 
  geom_histogram(aes(y = ..density..), binwidth = 5) + 
  geom_density()
Тьерри
источник
28

Вот как я это делаю:

foo <- rnorm(100, mean=1, sd=2)
hist(foo, prob=TRUE)
curve(dnorm(x, mean=mean(foo), sd=sd(foo)), add=TRUE)

Бонусное упражнение - сделать это с помощью пакета ggplot2 ...

Джон Джонсон
источник
Однако, если вам нужно что-то смещенное, вы можете либо выполнить приведенный выше пример плотности, преобразовать данные (например, foo.log & lt; - log (foo) и попробовать то, что указано выше), либо попробовать установить искаженное распределение, например гамма или логнормальная (логнормальная эквивалентна взятию журнала и подгонке нормального, кстати).
Джон Джонсон,
2
Но для этого по-прежнему необходимо сначала оценить параметры вашего распределения.
Дирк Эддельбюттель,
Это немного уходит от простого обсуждения R, поскольку мы углубляемся в теоретическую статистику, но вы можете попробовать эту ссылку для Gamma: en.wikipedia.org/wiki/Gamma_distribution#Parameter_estimation Для логнормального просто возьмите журнал (предполагая все данные положительные) и работают с данными, преобразованными в журнал. Для чего-нибудь более интересного, я думаю, вам придется поработать с учебником по статистике.
Джон Джонсон,
3
Я думаю, вы неправильно понимаете, как и исходный плакат, и все другие ответы вполне довольны использованием непараметрических оценок - таких как гистограмма старой школы или несколько более современная оценка плотности на основе данных. Параметрические оценки хороши, если у вас есть веские основания подозревать наличие распределения. Но здесь все было не так.
Дирк Эддельбюттель,
11

Дирк объяснил, как построить функцию плотности на гистограмме. Но иногда вам может потребоваться более сильное предположение об искаженном нормальном распределении и построить график вместо плотности. Вы можете оценить параметры распределения и построить его с помощью пакета sn :

> sn.mle(y=c(rep(65, times=5), rep(25, times=5), rep(35, times=10), rep(45, times=4)))
$call
sn.mle(y = c(rep(65, times = 5), rep(25, times = 5), rep(35, 
    times = 10), rep(45, times = 4)))

$cp
    mean     s.d. skewness 
41.46228 12.47892  0.99527

График распределенных данных с перекосом нормалей

Вероятно, это лучше работает с данными, которые более искажены:

Еще один искаженный сюжет

fmark
источник
3

У меня была та же проблема, но решение Дирка, похоже, не сработало. Я получал это предупреждение каждый раз

"prob" is not a graphical parameter

Я прочитал ?histи нашел околоfreq: a logical vector set TRUE by default.

код, который работал у меня,

hist(x,freq=FALSE)
lines(density(x),na.rm=TRUE)
Матиас Андина
источник