已知条件概率,反推设计值
已知: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_greaterv <- 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_greatery <- 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")