编辑代码

png(file = "/tmp/output.png")

# Load the data
data(faithful)

# Create the histogram for the eruptions data
hist_data <- hist(faithful$eruptions, breaks=seq(min(faithful$eruptions), max(faithful$eruptions), length.out=30), plot=FALSE)

# Calculate density and scale it to match the histogram frequency scale
density_data <- density(faithful$eruptions)
density_data$y <- density_data$y * max(hist_data$counts) / max(density_data$y)

# Plot the histogram using barplot for better control
barplot(hist_data$counts, names.arg=hist_data$mids, col="grey", border="black", xlab="Duration", ylab="Frequency", main="")

# Add density line
lines(density_data, col="blue", lwd=2)

# Add normal distribution line
x <- seq(min(faithful$eruptions), max(faithful$eruptions), length.out=length(density_data$x))
y <- dnorm(x, mean=mean(faithful$eruptions), sd=sd(faithful$eruptions))
y <- y * max(hist_data$counts) / max(y)
lines(x, y, col="red", lwd=2)

# Add jittered points over histogram bars
jitter_width <- (max(faithful$eruptions) - min(faithful$eruptions)) / length(hist_data$counts)
for(i in 1:length(hist_data$counts)) {
  points(jitter(rep(hist_data$mids[i], hist_data$counts[i]), amount=jitter_width), rep(i - 0.5, hist_data$counts[i]), pch=16, col="black")
}

# Add the p-value text
p_value <- 2.32119725434119e-14
text(max(hist_data$mids), max(hist_data$counts), labels=bquote("P = " * .(format(p_value, scientific=TRUE))), pos=4)

dev.off()