Ajustando uma curva de densidade a um histograma em R

91

Existe uma função em R que ajusta uma curva a um histograma?

Digamos que você tenha o seguinte histograma

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

Parece normal, mas está distorcido. Eu quero ajustar uma curva normal que é enviesada para envolver este histograma.

Esta pergunta é bastante básica, mas não consigo encontrar a resposta para R na internet.

r histogram curve-fitting r-faq user5243421
fonte
Você deseja encontrar m e s de forma que a distribuição gaussiana N (m, s) se ajuste aos seus dados?
SteinNorheim
Não tenho certeza do que isso significa ...> _>
user5243421
10
@mathee: Eu acho que ele quer dizer m = média e s = desvio padrão. Distribuição gaussiana é outro nome para distribuição normal.
Peter Mortensen

Respostas:

154

Se entendi sua pergunta corretamente, então você provavelmente quer uma estimativa de densidade junto com o histograma:

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

Edite um bom tempo depois:

Aqui está uma versão um pouco mais arrumada:

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)

junto com o gráfico que produz:

insira a descrição da imagem aqui

Dirk Eddelbuettel
fonte
3
+1 - você também pode fazer ao contrário, ou seja, ajustar o gráfico de densidade para ajustar o histograma?
vonjd
2
Eu sugiro dar parâmetro adicional para lines(density(X,na.rm= TRUE)o vetor pode conter valores NA.
Anirudh
30

Isso é fácil com 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()

ou para imitar o resultado da solução de Dirk

ggplot(dataset, aes(x = X)) + 
  geom_histogram(aes(y = ..density..), binwidth = 5) + 
  geom_density()
Thierry
fonte
28

É assim que eu faço:

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

Um exercício bônus é fazer isso com o pacote ggplot2 ...

John johnson
fonte
No entanto, se quiser algo distorcido, você pode fazer o exemplo de densidade acima, transformar seus dados (por exemplo, foo.log & lt; - log (foo) e tentar o acima) ou tentar ajustar uma distribuição distorcida, como o gama ou lognormal (lognormal é equivalente a tomar o log e ajustar um normal, aliás).
John Johnson
2
Mas isso ainda requer estimar os parâmetros de sua distribuição primeiro.
Dirk Eddelbuettel
Isso vai um pouco longe de simplesmente discutir R, já que estamos nos aprofundando em estatísticas teóricas, mas você pode tentar este link para o Gama: en.wikipedia.org/wiki/Gamma_distribution#Parameter_estimation Para lognormal, apenas pegue o log (assumindo todos os dados são positivos) e trabalhar com dados transformados em log. Para qualquer coisa mais sofisticada, acho que você teria que trabalhar com um livro de estatística.
John Johnson
3
Acho que você não entendeu como tanto o autor da postagem original quanto todas as outras respostas se contentam em usar estimativas não paramétricas - como um histograma da velha escola ou uma estimativa de densidade baseada em dados um pouco mais moderna. As estimativas paramétricas são ótimas se você tiver bons motivos para suspeitar de uma distribuição. Mas não foi esse o caso aqui.
Dirk Eddelbuettel
11

Dirk explicou como representar graficamente a função de densidade no histograma. Mas às vezes você pode querer ir com a suposição mais forte de uma distribuição normal distorcida e plotar isso em vez da densidade. Você pode estimar os parâmetros da distribuição e plotá-los usando o pacote 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

Gráfico de dados distribuídos com inclinação normal

Isso provavelmente funciona melhor com dados que são mais normais:

Outro gráfico enviesado normal

fmark
fonte
3

Eu tive o mesmo problema, mas a solução de Dirk não parecia funcionar. Eu recebia essa mensagem de aviso todas as vezes

"prob" is not a graphical parameter

Eu li ?histe descobri sobrefreq: a logical vector set TRUE by default.

o código que funcionou para mim é

hist(x,freq=FALSE)
lines(density(x),na.rm=TRUE)
Matias Andina
fonte