Dopasowanie krzywej gęstości do histogramu w R.

92

Czy w R jest funkcja dopasowująca krzywą do histogramu?

Powiedzmy, że masz następujący histogram

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

Wygląda normalnie, ale jest przekrzywiony. Chcę dopasować normalną krzywą, która jest skośna, aby zawijać się wokół tego histogramu.

To pytanie jest raczej podstawowe, ale nie mogę znaleźć odpowiedzi na R w Internecie.

r histogram curve-fitting r-faq

— user5243421
źródło

Czy chcesz znaleźć mi s takie, że rozkład Gaussa N (m, s) pasuje do twoich danych?

— SteinNorheim

Nie jestem pewien, co to oznacza ...> _>

— user5243421

10

@mathee: Myślę, że ma on na myśli m = średnia, a s = odchylenie standardowe. Rozkład Gaussa to inna nazwa dystrybucji normalnej.

— Peter Mortensen

155

Jeśli dobrze rozumiem twoje pytanie, prawdopodobnie chcesz oszacować gęstość wraz z histogramem:

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

Edytuj długo później:

Oto nieco bardziej wystrojona wersja:

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)

wraz z generowanym wykresem:

wprowadź opis obrazu tutaj

— Dirk Eddelbuettel
źródło

3

+1 - czy można to też zrobić na odwrót, tj. Dopasować wykres gęstości do histogramu?

— vonjd

2

Sugeruję podanie dodatkowego parametru, lines(density(X,na.rm= TRUE)ponieważ wektor może zawierać wartości NA.

— Anirudh

30

Z ggplot2 jest to łatwe

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()

lub aby naśladować wynik z rozwiązania Dirka

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

— Thierry
źródło

28

Oto jak to robię:

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

Dodatkowym ćwiczeniem jest zrobienie tego z pakietem ggplot2 ...

— John Johnson
źródło

Jeśli jednak chcesz czegoś, co jest wypaczone, możesz albo wykonać przykład gęstości z góry, przekształcić dane (np. Foo.log & lt; - log (foo) i wypróbować powyższe) lub spróbować dopasować skośną dystrybucję, na przykład gamma lub lognormal (lognormal jest równoznaczne z pobraniem logu i dopasowaniem normalnego, btw).

— John Johnson,

2

Ale to nadal wymaga najpierw oszacowania parametrów twojej dystrybucji.

— Dirk Eddelbuettel

To jest trochę dalekie od zwykłego omawiania R, ponieważ zagłębiamy się bardziej w statystyki teoretyczne, ale możesz spróbować tego linku do Gamma: en.wikipedia.org/wiki/Gamma_distribution#Parameter_estimation Dla lognormal, po prostu weź dziennik (zakładając wszystkie dane są dodatnie) i pracować z danymi przekształconymi w dziennik. Wydaje mi się, że dla kogoś bardziej wyszukanego musiałbyś popracować z podręcznikiem statystyki.

— John Johnson,

3

Myślę, że źle rozumiesz, że zarówno oryginalny plakat, jak i wszystkie inne odpowiedzi są całkiem zadowolone z użycia szacunków nieparametrycznych - takich jak oldschoolowy histogram lub nieco bardziej nowoczesne, zagęszczone oszacowanie oparte na danych. Szacunki parametryczne są świetne, jeśli masz dobry powód, by podejrzewać rozkład. Ale tak nie było w tym przypadku.

— Dirk Eddelbuettel

11

Dirk wyjaśnił, jak wykreślić funkcję gęstości na histogramie. Ale czasami możesz chcieć pójść z silniejszym założeniem skośnego rozkładu normalnego i wykreślić to zamiast gęstości. Możesz oszacować parametry dystrybucji i wykreślić ją za pomocą pakietu 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

Wykres danych z rozkładem normalnym i skośnym

Prawdopodobnie działa to lepiej w przypadku danych, które są bardziej normalne:

Kolejna normalna fabuła

— fmark
źródło

3

Miałem ten sam problem, ale rozwiązanie Dirka wydawało się nie działać. Za każdym razem otrzymywałem to ostrzeżenie

"prob" is not a graphical parameter

Przeczytałem ?histi znalazłemfreq: a logical vector set TRUE by default.

kod, który działał dla mnie, to

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

— Matias Andina
źródło