基于R语言的高斯拟合可视化扩展及ggplot2适配技术问询

Hey there! Let's work through your two requirements to polish up your Gaussian fit visualization. I'll break this into replacing the base R plot with ggplot2, extending the curve beyond day 35, and adding uncertainty bands to show time-dimension variability.

Step 1: Prep Data & Refine the Fitting Function #

First, we'll adjust your existing fit function to return the hessian matrix—this lets us calculate parameter uncertainty, which is key for the error bands later. We'll also load the required packages:

# Install required packages if you haven't already
install.packages(c("ggplot2", "MASS"))

# Load packages
library(ggplot2)
library(MASS)

# Your original data
x <- c(1:35)
y <- c(221,88,76,203,233,228,288,498,428,443,570,640,1145,1326,1598, 529,2076,2249,2116,2795,2853,2470,2989,2648,4480,4670,4821, 3957,3780,3612,3491,4492,4401,3651,3815)
data <- data.frame(x, y)

# Updated fit function with hessian output for uncertainty calculations
fitG <- function(x, y, mu, sig, scale) {
  f <- function(p) {
    d <- p[3] * dnorm(x, mean = p[1], sd = p[2])
    sum((d - y)^2)
  }
  # Request hessian matrix to compute parameter standard errors
  optim(c(mu, sig, scale), f, hessian = TRUE)
}

# Run the fit with your initial guesses
fitP <- fitG(data$x, data$y, 35, 1, 6000)

Step 2: Generate Extended Fit Data & Uncertainty Bands #

Next, we'll create an extended time range (I'll go up to day 45, but you can adjust this), calculate the fitted Gaussian curve for these days, and simulate uncertainty bands using the parameter covariance from the hessian:

# Extract fitted parameters and their covariance matrix
params <- fitP$par
param_cov <- solve(fitP$hessian)  # Covariance matrix of parameters

# Create extended time sequence (adjust the upper limit as needed)
x_extended <- seq(1, 45, by = 0.5)

# Helper function to calculate Gaussian fit values
gaussian_fit <- function(x, mu, sig, scale) {
  scale * dnorm(x, mean = mu, sd = sig)
}

# Calculate base fitted values for the extended range
fitted_vals <- gaussian_fit(x_extended, params[1], params[2], params[3])

# Simulate 1000 sets of parameters from their multivariate normal distribution
set.seed(123)  # Ensure reproducible results
sim_params <- mvrnorm(n = 1000, mu = params, Sigma = param_cov)

# Generate fitted curves for each simulated parameter set
sim_fits <- apply(sim_params, 1, function(p) gaussian_fit(x_extended, p[1], p[2], p[3]))

# Calculate 95% confidence bands (2.5th and 97.5th percentiles)
lower_band <- apply(sim_fits, 1, quantile, 0.025)
upper_band <- apply(sim_fits, 1, quantile, 0.975)

# Combine all fit data into a single data frame for ggplot
fit_data <- data.frame(
  x = x_extended,
  fitted = fitted_vals,
  lower = lower_band,
  upper = upper_band
)

Step 3: Plot with ggplot2 #

Now we'll build the ggplot with your original data points, the extended fitted curve, and the uncertainty bands. We'll add custom labels and a clean theme for readability:

ggplot() +
  # Plot original data points
  geom_point(data = data, aes(x = x, y = y), color = "darkslateblue", size = 2) +
  # Plot the main fitted Gaussian curve
  geom_line(data = fit_data, aes(x = x, y = fitted), color = "firebrick", linewidth = 1.2) +
  # Add 95% uncertainty bands (semi-transparent ribbon)
  geom_ribbon(data = fit_data, aes(x = x, ymin = lower, ymax = upper), 
              fill = "firebrick", alpha = 0.2) +
  # Customize labels and theme
  labs(
    title = "Gaussian Fit of COVID Case Data",
    x = "Day",
    y = "Reported Cases",
    caption = "Data source: Italian Civil Protection"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    plot.caption = element_text(hjust = 0, size = 10, color = "gray50")
  )

Key Notes: #

• Extending the curve: Change the upper value in x_extended (e.g., seq(1, 50, by = 1)) to cover more days beyond 35.
• Uncertainty bands: The 95% ribbon shows the range of plausible fit curves given the uncertainty in your Gaussian parameters (mean, SD, scale). This directly translates to uncertainty in the time dimension (e.g., when the peak might occur, how the curve will decay after day 35).
• Reproducibility: The set.seed(123) ensures your simulated bands are the same every time you run the code.

内容的提问来源于stack exchange，提问作者Raffaello Nardin