已知条件概率,反推设计值

已知:P(Y>y|X>x), P(X>x)=t的情况下,求P(Y>y|X>x)为某一概率时Y对应的设计值。

###已知条件概率求边缘设计值#####


# 定义分布参数
# PearsonⅢ分布参数
shape_p3 <- 0.2557214
location_p3 <- 12.87857
scale_p3 <- 67.14898

# gamma分布参数
shape_gamma <- 1.7192805823
rate_gamma <- 0.0001989357

# Frank copula参数
theta_frank <- -1.365

# 创建Frank copula对象
frank_cop <- frankCopula(param = theta_frank)

# 定义P(X > x)的值
prob_x_greater_values <- c(0.1, 0.05, 0.02)

# 定义P(Y > y|X > x)的值
conditional_prob_values <- c(0.5, 0.4, 0.3, 0.2, 0.1)

# 定义求解y的函数
solve_y <- function(target_joint_prob, prob_x_greater, cop) {
  f <- function(y) {
    # 计算 P(Y > y)
    prob_y_greater <- 1 - pgamma(y, shape = shape_gamma, rate = rate_gamma)
    u <- 1-prob_x_greater
    v <- 1-prob_y_greater
    # 计算联合概率分布 P(X > x, Y > y)
    joint_prob <-1-u-v+ pCopula(c(u, v), cop)
    return(joint_prob - target_joint_prob)
  }
  
  intervals <- list(c(0, 1e6))
  for (interval in intervals) {
    result <- tryCatch({
      uniroot(f, interval = interval)
    }, error = function(e) {
      NULL
    })
    if (!is.null(result)) {
      return(result$root)
    }
  }
  cat("无法找到合适的y值,目标联合概率为", target_joint_prob, "\n")
  return(NA)
}
result_matrix <- matrix(NA, nrow = length(prob_x_greater_values), ncol = length(conditional_prob_values))
# 遍历不同的P(X > x)和P(Y > y|X > x)组合
for (i in seq_along(prob_x_greater_values)) {
  prob_x_greater <- prob_x_greater_values[i]
  for (j in seq_along(conditional_prob_values)) {
    conditional_prob <- conditional_prob_values[j]
    # 计算P(Y > y, X > x)
    target_joint_prob <- conditional_prob * prob_x_greater
    y <- solve_y(target_joint_prob, prob_x_greater, frank_cop)
    result_matrix[i, j] <- y
  }
}

# 保存结果到.xlsx文件
file_path <- file.path(getwd(), "result_table.xlsx")
write.xlsx(result_matrix, file = file_path)
cat("结果已保存到", file_path, "\n")

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值