已知: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")