diff --git a/R/BaconMethods.R b/R/BaconMethods.R
index dacfc70..7c4a608 100644
--- a/R/BaconMethods.R
+++ b/R/BaconMethods.R
@@ -76,19 +76,26 @@ setMethod("traces", "Bacon", function(object, burnin=TRUE, index=1){
     else
         gstraces <- object@traces[-c(1:object@nburnin),,index]
 
-    thetahat <- estimates(object)[index,]
-
-    op <- par(mfcol=c(3, 3), mar=c(2,4,2,2))
-    for(i in 1:9) {
-        if(burnin)
-            plot(1:object@niter, gstraces[,i],
-                 ylab=colnames(gstraces)[i], type="l", xlab="", main = "", lwd=0.3)
-        else
-            plot((object@nburnin+1):object@niter, gstraces[,i],
-                 ylab=colnames(gstraces)[i], type="l", xlab="", main = "", lwd=0.3)
-        abline(h=thetahat[i], col=2, lwd=2)
-    }
-    par(op)
+    thetahat <- as.data.frame(estimates(object)[index,]) %>% 
+        rownames_to_column("variable") %>%
+        rename(value = "estimates(object)[index, ]")
+    g <- ggplot(gstraces_melted, aes(x = iteration, y = value)) +
+    geom_line() +
+    geom_hline(data = thetahat_df, aes(yintercept=value), color = "red") +
+    facet_wrap(~variable, scales = "free_y", strip.position = "left") +
+    scale_x_continuous(labels = c("",1000, "", 3000, "", 5000)) +
+    theme_cowplot(font_size = 12) +
+    theme(strip.background = element_blank(),
+          strip.placement = "outside",
+          plot.margin = margin(0.15, 0.25, 0.15, 0.15, "in")) +
+    xlab("Iteration") +
+    ylab("Trace")
+  
+    if (!burnin) {
+    g <- g + xlim(c((object@nburnin+1), object@niter))
+    } 
+  
+    return(g)
 })
 
 ##' @rdname posteriors-methods
@@ -103,17 +110,26 @@ setMethod("posteriors", "Bacon", function(object, thetas, index, alphas, xlab, y
     if(xlab=="") xlab <- thetas[1]
     if(ylab=="") ylab <- thetas[2]
 
-    plot(gstraces, pch=20, xlab=xlab, ylab=ylab, bty='n',
-         main=c("median at:", round(estimates(object)[thetas], 3)))
-    points(estimates(object)[index, thetas], col=3, pch=17, cex=2)
-
-    for(alpha in alphas)
-        lines(ellipse(cov(gstraces), centre=colMeans(gstraces), level=alpha), col="blue", ...)
+    df <- data.frame(x = gstraces[,1], y = gstraces[,2])
+    est_df <- data.frame(x = estimates(object)[index, thetas[1]], y = estimates(object)[index, thetas[2]])
+    
+    # Plot using ggplot
+    p <- ggplot(df, aes(x = x, y = y)) +
+        geom_point(shape = 20) +
+        stat_ellipse(level = 0.95, col = "blue") +
+        stat_ellipse(level = 0.9, col = "blue") +
+        stat_ellipse(level = 0.75, col = "blue") +
+        labs(x = xlab,
+             y = ylab) + 
+        ggtitle(paste("median at:", round(estimates(object)[thetas], 3))) +
+        geom_point(data = est_df, aes(x = x, y = y), color = "red", shape = 17, size = 4)
+    
+    return(p)
 })
 
 ##' @rdname fit-methods
 ##' @aliases fit
-setMethod("fit", "Bacon", function(object, index, col="grey75", border="grey75", ...){
+setMethod("fit", "Bacon", function(object, index, ...){
     plotnormmix(tstat(object, corrected=FALSE)[, index], estimates(object)[index, ], ...)
 })
 
diff --git a/R/normmixture.R b/R/normmixture.R
index 9160aa5..90df7e8 100644
--- a/R/normmixture.R
+++ b/R/normmixture.R
@@ -74,12 +74,26 @@ dnormmix <- function(x, theta){
 plotnormmix <- function(x, theta, ...) {
     if(length(theta) %% 3 != 0)
         stop("Length of theta should be a multiple of three!")
-    hist(x=x, freq=FALSE, ...)
-    f <- function(x) dnormmix(x, theta)
-    curve(expr=f, add=TRUE, col=1, lwd=2)
-    ncomp <- length(theta)/3
-    for(k in 1:ncomp) {
-        f <- function(x) theta[k]*dnorm(x, mean=theta[k+3], sd=theta[k+6])
-        curve(expr=f, add=TRUE, col=k+1, lwd=2)
-    }
+    x <- data.frame(x = x)
+    theta <- data.frame(y = theta)
+    fit <- ggplot(x, aes(x=x , y = after_stat(density))) +
+    geom_histogram(fill = "grey", color="black", binwidth = 1.5) +
+    geom_line(aes(x=x, y =dnorm(x, mean(x), sd(x))), lwd=1) +
+    geom_line(aes(x=x, 
+              y=theta["p.0",]*dnorm(x, theta["mu.0",], theta["sigma.0",])), 
+              color="red",
+              lwd = 1.5) +
+    geom_line(aes(x=x, 
+              y=theta["p.1",]*dnorm(x, theta["mu.1",], theta["sigma.1",])), 
+              color="green",
+              lwd=1.5) +
+    geom_line(aes(x=x, 
+              y=theta["p.2",]*dnorm(x, theta["mu.2",], theta["sigma.2",])), 
+              color="blue",
+              lwd=1.5) +
+    theme_cowplot(font_size = 12) +
+    xlab("Test Statistics") +
+    ylab("Density")
+  
+    return(fit)
 }