Created
September 6, 2016 11:08
-
-
Save claczny/556f768d4c972cf4d62abcea2850b33c to your computer and use it in GitHub Desktop.
Fit and plot the fitted components of a multi-modal distribution assuming normal distributions as the components.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(ggplot2) | |
library(plyr) | |
library(mixtools) | |
### | |
# GLOBAL THEME AND GLOBAL AESTHETICS | |
### | |
old <- theme_set(theme_bw() + | |
theme(text = element_text(size=12), | |
axis.title = element_text(size = 14, face="bold"), | |
title = element_text(face = "bold") | |
) | |
) | |
update_geom_defaults("point", list(size = 2.5)) | |
### | |
# FUNCTIONS ### | |
### | |
# Motivated by http://tinyheero.github.io/2016/01/03/gmm-em.html | |
#' Plot a Mixture Component | |
#' | |
#' @param x Input data. | |
#' @param mu Mean of component. | |
#' @param sigma Standard of component. | |
#' @param lam Mixture weight of component. | |
plot_mix_comps <- function(x, mu, sigma, lam) { | |
lam * dnorm(x, mu, sigma) | |
} | |
### | |
# DATA ##### | |
### | |
# data(faithful) | |
# attach(faithful) | |
df <- waiting | |
# Fit the normal mixture(s) | |
mixmdl <- normalmixEM(df) | |
# Base graphics | |
# plot(mixmdl,which=2, breaks=200) | |
# lines(density(df), lty=2, lwd=2) | |
# Using ggplot | |
# Motivated by http://tinyheero.github.io/2016/01/03/gmm-em.html | |
p <- data.frame(x = df) %>% | |
ggplot(aes(x=x)) + | |
geom_histogram(aes(x, ..density..), binwidth = 1, colour = "black", | |
fill = "white") + | |
stat_function(fun = plot_mix_comps, aes(colour="1"), | |
args = list(mixmdl$mu[1], mixmdl$sigma[1], | |
lam = mixmdl$lambda[1]), lwd = 1.5) + | |
stat_function(fun = plot_mix_comps, aes(colour="2"), | |
args = list(mixmdl$mu[2], mixmdl$sigma[2], | |
lam = mixmdl$lambda[2]), lwd = 1.5) + | |
geom_density(aes(x, ..density.., colour="joint"), size=1.5) + | |
scale_colour_manual("Component", values = c("1"="red", "2"="blue", "joint" = "black")) + | |
ylab("Density") + | |
xlab("Values") | |
print(p) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
DESRCIPTION
mixtools
provide plotting functionality for the fits, e.g., viaplot.mixEM
usingbase
graphics (AFAIK).This code is meant to demonstrate the use of
ggplot
to create a plot of the data histogram, the "joint" (prior) distribution, and the (posterior) distributions of the individual components, here two components.Should your data have more than two modes, you should add
stat_function()
entries according to the number of fitted components.