R语言实现均匀分布采样下多项式求解与3D密度图绘制

核心逻辑说明 #

• 约束化简：给定约束a + b < 1 - a + b可直接消去两侧b项，等价于a < 0.5，采样时无需重复判断不等式，直接按范围采样即可
• 方程化简：原方程展开后x的三次项会完全抵消，实际为关于x的一元二次方程，直接用求根公式计算即可，无需求助第三方多项式包，数值稳定性更高
• 结果存储：循环采样时动态存储所有有效实根，支持自定义theta值、采样次数t
• 3D绘图：采用三维核密度估计绘制交互式密度图，支持鼠标拖拽旋转查看分布

完整可运行代码 #

# 加载绘图依赖包，未安装请先运行 install.packages(c("rgl", "MASS"))
library(rgl)
library(MASS)

poly_solve_3d <- function(theta = 1, t = 1000, a_min = 0, b_min = 0) {
  # 预分配结果存储空间，提升大样本下的运行速度
  res <- data.frame(
    a = numeric(2*t),
    b = numeric(2*t),
    x = numeric(2*t)
  )
  valid_cnt <- 0

  for (i in 1:t) {
    # 采样满足约束的参数a,b：a<0.5，b默认取0到1-a范围保证s=1-a-b非负
    a <- runif(1, min = a_min, max = 0.5)
    b <- runif(1, min = b_min, max = 1 - a)
    s <- 1 - a - b

    # 计算一元二次方程 A*x² + B*x + C = 0 的系数
    A <- a - b - theta * s
    B <- theta * s^2 - 2*a*s
    C <- a * s^2

    roots <- c()
    # 处理退化为一次方程的边界情况
    if (abs(A) < 1e-9) {
      if (abs(B) > 1e-9) roots <- -C/B
    } else {
      # 二次方程求根，考虑浮点误差修正判别式
      delta <- B^2 - 4*A*C
      delta <- max(delta, 0)
      x1 <- (-B + sqrt(delta))/(2*A)
      x2 <- (-B - sqrt(delta))/(2*A)
      roots <- c(x1, x2)
    }

    # 过滤无效值，存入结果
    roots <- roots[is.finite(roots)]
    if (length(roots) > 0) {
      root_num <- length(roots)
      res[(valid_cnt+1):(valid_cnt+root_num), ] <- data.frame(
        a = rep(a, root_num),
        b = rep(b, root_num),
        x = roots
      )
      valid_cnt <- valid_cnt + root_num
    }
  }

  # 裁剪空行
  res <- res[1:valid_cnt, ]

  # 绘制3D交互式密度图
  if (valid_cnt > 20) {
    kde_est <- kde3d(res$a, res$b, res$x, n = 50)
    persp3d(
      kde_est$x, kde_est$y, kde_est$z,
      col = "lightblue", alpha = 0.6,
      xlab = "a", ylab = "b", zlab = "x",
      main = paste0("3D密度分布 | theta = ", theta)
    )
    # 若需要叠加原始散点可取消下一行注释
    # points3d(res$a, res$b, res$x, col = rgb(1,0,0,0.2), size = 2)
  }

  return(res)
}

# 调用示例：theta设为1.5，采样2000次
result <- poly_solve_3d(theta = 1.5, t = 2000)

自定义调整说明 #

• 若仅需保留正/负x解，在roots计算完成后增加子集筛选即可：正解添加roots <- roots[roots > 1e-6]，负解添加roots <- roots[roots < -1e-6]
• 若需要调整a、b的采样范围，直接修改函数的a_min/b_min参数，或修改b采样的上界逻辑即可
• 大样本量（t>10000）场景下可直接调整预分配空间的大小，避免内存动态扩容拖慢速度
• 生成的rgl 3D图像支持鼠标拖拽旋转、滚轮缩放，可调整alpha参数修改曲面透明度查看内部密度分布

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