如何获取R语言中density()生成的核密度图拐点的x值？

Got it, let's walk through how to pull the x-values of inflection points from a kernel density curve generated by R's density() function. Since you already know how to find extrema, this builds on similar numeric reasoning—inflection points are where the second derivative of the density curve equals 0 (and the concavity flips), so we just need to extend that logic to second derivatives.

Here's a step-by-step, reproducible approach:

1. Generate the Kernel Density Estimate #

First, create your density object using your data. We'll use mtcars$mpg as an example:

# Load sample data and generate density curve
data <- mtcars$mpg
dens <- density(data, n = 1024)  # Use a higher n for more precise calculations

Note: Increasing n (the number of x-points generated) will improve the accuracy of our derivative calculations.

2. Calculate the Second Derivative Numerically #

Since density() returns a discrete set of x and y (density) values, we'll use numeric differentiation to find the second derivative:

# Extract the evenly spaced x-values and density values
x_vals <- dens$x
y_vals <- dens$y

# Calculate the step size between x-points (uniform, so dx is constant)
dx <- diff(x_vals)[1]

# First derivative: difference in y divided by difference in x
first_deriv <- diff(y_vals) / dx
# The first derivative corresponds to midpoints between original x-values
first_deriv_x <- x_vals[-1] - dx/2

# Second derivative: difference in first derivative divided by dx
second_deriv <- diff(first_deriv) / dx
# The second derivative corresponds to midpoints of the first derivative x-values
second_deriv_x <- first_deriv_x[-1] - dx/2

3. Locate Inflection Points (Where Second Derivative Flips Sign) #

Inflection points occur where the second derivative changes sign (from positive to negative, or vice versa). We'll use linear interpolation to find the exact x-value where the second derivative crosses 0:

# Find indices where the sign of the second derivative changes
sign_changes <- which(diff(sign(second_deriv)) != 0)

# Interpolate to find the exact x-values where second derivative = 0
inflection_x <- sapply(sign_changes, function(i) {
  # Grab the endpoints of the interval where the sign changes
  x1 <- second_deriv_x[i]
  x2 <- second_deriv_x[i+1]
  y1 <- second_deriv[i]
  y2 <- second_deriv[i+1]
  
  # Linear interpolation to solve for x when y=0
  x1 + (0 - y1) * (x2 - x1) / (y2 - y1)
})

4. Verify with a Plot #

To make sure we got it right, plot the density curve and overlay the inflection points:

plot(dens, main = "Kernel Density of mpg with Inflection Points")
abline(v = inflection_x, col = "red", lty = 2, lwd = 1.5)
legend("topright", legend = "Inflection Points", col = "red", lty = 2)

Key Notes #

• Bandwidth Impact: Just like with extrema, the kernel bandwidth (bw parameter in density()) will change the number and position of inflection points. Wider bandwidths create smoother curves with fewer inflection points; narrower bandwidths create more wiggly curves with more inflection points.
• Alternative with numDeriv Package: For a more streamlined (but package-dependent) approach, use the numDeriv package to compute second derivatives directly:

  library(numDeriv)
  
  # Create an interpolated function for the density curve
  dens_fun <- approxfun(dens$x, dens$y)
  
  # Define a function to return the second derivative at a given x
  second_deriv_fun <- function(x) hessian(dens_fun, x)
  
  # Evaluate second derivative across a dense grid of x-values
  x_grid <- seq(min(dens$x), max(dens$x), length.out = 2000)
  second_deriv_vals <- sapply(x_grid, second_deriv_fun)
  
  # Find sign changes and interpolate inflection points (same as step 3)
  sign_changes <- which(diff(sign(second_deriv_vals)) != 0)
  inflection_x <- sapply(sign_changes, function(i) {
    x1 <- x_grid[i]
    x2 <- x_grid[i+1]
    y1 <- second_deriv_vals[i]
    y2 <- second_deriv_vals[i+1]
    x1 + (0 - y1) * (x2 - x1) / (y2 - y1)
  })
  

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