如何在R中绘制Birnbaum-Saunders分布的风险函数？

你提到的分母构造问题可以直接用R内置函数解决：pnorm()就是R原生的标准正态分布累积分布函数，对应的概率密度函数为dnorm()，直接调用即可完成公式构造。

实现步骤 #

• 第一步：自定义Birnbaum-Saunders（BS）分布的风险函数

# 自定义BS分布风险函数，默认尺度参数β=1
h_bs <- function(t, alpha) {
  # 计算参数ε(t)
  epsilon <- sqrt(t) - 1/sqrt(t)
  # 计算BS分布的概率密度f(t)
  ft <- dnorm(epsilon/alpha) * (sqrt(t) + 1/sqrt(t)) / (2 * alpha * t)
  # 计算BS分布的生存函数S(t)=1-F(t)，即你需要的分母部分
  st <- 1 - pnorm(epsilon/alpha)
  # 返回风险函数计算结果
  return(ft/st)
}

• 第二步：生成绘图用的自变量序列，注意t>0，需避开t=0的奇点

# 生成t的取值序列，范围0.01到10，步长0.01，可根据需求调整范围
t_seq <- seq(0.01, 10, by = 0.01)

• 第三步：设定参数α的取值，可设置多个值对比不同参数下的曲线差异

# 设定多组α参数测试
alpha_list <- c(0.5, 1, 2)
# 批量计算对应风险函数值
h_values <- sapply(alpha_list, function(a) h_bs(t_seq, a))

• 第四步：调用plot()函数绘制图形

# 绘制第一条α对应的风险曲线
plot(t_seq, h_values[,1], type = "l", lwd = 2, col = "red",
     xlab = "t", ylab = "h(t;α)", main = "Birnbaum-Saunders分布风险函数")
# 追加其他α的曲线
lines(t_seq, h_values[,2], lwd = 2, col = "blue")
lines(t_seq, h_values[,3], lwd = 2, col = "green")
# 添加图例
legend("topright", legend = paste0("α=", alpha_list),
       col = c("red", "blue", "green"), lwd = 2)

注意事项 #

• 如果计算时出现NaN，可适当调大t的最小取值，避开奇点区域即可
• 如果你要绘制自己研究的类BS分布，只需修改自定义函数中ft和st的计算逻辑即可，核心都是用pnorm()调用标准正态CDF构造分母

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