File size: 3,168 Bytes
998922f | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 | library(ggplot2)
library(dplyr)
# Create a bimodal distribution function (mixture of two normals)
bimodal <- function(x) {
0.6 * dnorm(x, mean = -1, sd = 0.8) +
0.4 * dnorm(x, mean = 2, sd = 1.2)
}
# Generate data points
x <- seq(-4, 4, length.out = 1000)
y <- bimodal(x)
df <- data.frame(x = x, y = y)
# Function to numerically find quantiles for our bimodal distribution
find_quantile <- function(p) {
# Create a function that we want to minimize
f <- function(x) {
integrate(bimodal, -Inf, x)$value - p
}
# Use uniroot to find where this function equals 0
uniroot(f, interval = c(-10, 10))$root
}
# Calculate quantiles
probs <- seq(0.05, 0.95, length.out = 15)
quantiles <- c(-20, sapply(probs, find_quantile), 20)
# Create data frame for the clipped vertical lines
line_data <- data.frame()
for (q in quantiles[2:(length(quantiles) - 1)]) { # Skip the ±20 points
y_at_q <- bimodal(q)
line_data <- rbind(
line_data,
data.frame(
x = q,
y_start = 0,
y_end = y_at_q
)
)
}
# Create data for filled intervals
interval_data <- data.frame()
for (i in 1:(length(quantiles) - 1)) {
x_seq <- seq(max(-4, quantiles[i]),
min(4, quantiles[i + 1]),
length.out = 100
)
interval_data <- rbind(
interval_data,
data.frame(
x = x_seq,
y = bimodal(x_seq),
group = i
)
)
}
# Create the plot
p <- ggplot() +
# Add filled intervals
geom_ribbon(
data = interval_data,
aes(x = x, ymin = 0, ymax = y, group = group),
fill = "lightgray",
alpha = 0.3
) +
# Add the distribution curve
geom_line(
data = df, aes(x = x, y = y),
color = "black", size = 1
) +
# Add thin dashed vertical lines for quantiles
geom_segment(
data = line_data,
aes(
x = x, xend = x,
y = y_start, yend = y_end
),
color = "darkgray",
linetype = "dashed",
size = 0.3,
alpha = 0.7
) +
# Add quantile labels
geom_text(
data = data.frame(
x = quantiles[2:(length(quantiles) - 1)], # Skip the ±20 points
y = rep(-0.01, 15),
label = paste0("", 1:15)
),
aes(x = x, y = y, label = label),
angle = -30,
vjust = 1,
size = 4
) +
# Add infinity labels
geom_text(
data = data.frame(
x = c(-4, 4),
y = rep(-0.01, 2),
label = c("0", "16")
),
aes(x = x, y = y, label = label),
angle = -30,
vjust = 1,
size = 4
) +
# Customize the theme and labels
theme_bw(base_size = 16) +
theme(
axis.line = element_line(color = "black"),
plot.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank()
) +
labs(
x = "x",
y = "Density"
) +
# Set the axis limits
scale_x_continuous(limits = c(-4, 4)) +
ylim(-0.02, 0.32) +
# Add theme customizations
theme(
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.grid.minor = element_blank()
)
# Save the plot as PDF
ggsave(
"bimodal_quantiles.pdf",
plot = p,
width = 10, # Width in inches
height = 7, # Height in inches
device = "pdf",
dpi = 300 # High resolution
)
|