par(mfrow = c(2,5))
for(j in 1:6) {
#    hist(output[1:200,j], main = names(output)[j], freq=F)
      xmean = summary(Rfit)@coefs[j,1]
      xsd = summary(Rfit)@coefs[j,2]
      xmin = xmean -3*xsd; xmax = xmean + 3*xsd
      xvals = seq(xmin, xmax, (xmax-xmin)/100)
      yvals = dnorm(xvals, mean=xmean, sd=xsd)
    minx = min(density(output[,j])$x, density(output2[,j])$x, xvals)
    maxx = max(density(output[,j])$x, density(output2[,j])$x, xvals)
    miny = 0
    maxy = max(density(output[,j])$y, density(output2[,j])$y, yvals)
    plot(xvals, yvals, col = "black", type = "l", lwd = 2, main = names(output)[j], xlim = c(minx, maxx), ylim = c(miny, maxy), xlab = "estimate", ylab = "density")
    lines(density(output[,j]), col = "red", lwd = 2)
    lines(density(output2[,j]), col = "blue", lwd = 2)
}
for(j in 7:10) {
    minx = min(density(output[,j])$x, density(output2[,j])$x)
    maxx = max(density(output[,j])$x, density(output2[,j])$x)
    miny = 0
    maxy = max(density(output[,j])$y, density(output2[,j])$y)
    plot(density(output[,j]), col = "red", type = "l", lwd = 2, main = names(output)[j], xlim = c(minx, maxx), ylim = c(miny, maxy), xlab = "estimate", ylab = "density")
    lines(density(output2[,j]), col = "blue", lwd = 2)
}
legend(x="topleft", col = c("black", "blue", "red"), lwd = 2)