计算一个联合的条件概率p.d.f
Calculate the conditional probability of a joint p.d.f
我有关节p.d.f。
我现在正在比较条件概率的理论值和我的经验值
运行 Monte Carlo 方法。
我必须重复 10,000、100,000 和 1,000,000 次抽奖。请问如何在R代码中复制?
此外,对于最后一步,条件概率,
我正在使用 Monte Carlo 方法来做到这一点。
是否有任何 R 代码可用于多元制服
分布计算条件概率?
如有任何建议,我们将不胜感激!谢谢!!
我的代码如下:
# f(x.y) = (1/4)*xy, 0<x<2, 0<y<2
# Find P(A) = P(X>1)
f <- function(x){(1/2)*x} # Marginal P(X)
probE <-integrate(f, lower = 1, upper = 2)
cat('\n Pr[ 1 < X ] is \n')
print(probE)
n <- 10000
x<-runif(n, 1,2)
probE.MC <- ((2-1)/n)*sum((1/2)*x)
cat('\n Monte Carlo Pr[1< X ] =',probE.MC,'\n')
# Find P(B) = P(Y<1)
f <- function(y){(1/2)*y} # Marginal P(Y)
probB <-integrate(f, lower = 0, upper = 1)
cat('\n Pr[ 1< Y ] is \n')
probB
typeof(probB)
n <- 10000
y<-runif(n, 0,1)
probB.MC <- ((1-0)/n)*sum((1/2)*y)
cat('\n Monte Carlo Pr[Y < 1] =',probB.MC,'\n')
# Pr[A intersect B]
# P[X>1 and Y <1]
f <- function(x,y){return((1/4)*x*y)}
n <- 100000
a11<-1; a12 <-2; a21 <- 0; a22 <-1
x <-runif(n, a11, a12)
y <- runif(n,a21, a22)
probMC <- ((a12-a11)*(a22-a21)/n)*sum(f(x,y))
probMC
typeof(probMC)
# P[A|B] = p[A intersect B]/ P(B)
probAB <- probMC/probB
首先,我重新格式化了函数并给它们单独命名。
fX <- function(x) {
0.5 * x # Marginal P(X)
}
# Find P(B) = P(Y<1)
fY <- function(y) {
0.5 * y # Marginal P(Y)
}
fXY <- function(x, y) {
1 / length(x) * sum(0.25 * x * y) # joint X,Y
}
进行多次运行的最简单方法是将 MC 代码包装在一个 for 循环中,并将每个计算保存在一个数组中。然后,最后取存储值的平均值。
所以对于 P[A] 你有:
n <- 10000
probA.MC <- numeric(n) # create the array
for (i in 1:10000) {
x<-runif(n, 1,2)
probA.MC[i] <- ((2-1) / n) * sum(0.5 * x)
}
cat('\n Monte Carlo Pr[1 < X] =',mean(probA.MC),'\n')
(我假设 probE.MC 应该是 probA.MC。)结果是 Monte Carlo Pr[1 < X] = 0.7500088
。该代码类似于 P[B],结果为 Monte Carlo Pr[Y < 1] = 0.2499819
.
我们使用 fXY 作为联合概率。
n <- 10000
probMC <- numeric(n)
for (i in 1:10000) {
x <- runif(n, 1, 2)
y <- runif(n, 0, 1)
probMC[i] <- ((a12-a11) * (a22-a21)) * fXY(x, y)
}
cat('\n Monte Carlo Pr[X,Y] =',mean(probMC),'\n')
这个结果是Monte Carlo Pr[X,Y] = 0.1882728
。
您最后一次计算应该如下所示(注意积分结果中的 probB$ 值):
# P[A|B] = p[A intersect B]/ P(B)
probAB <- mean(probMC) / probB$value
print(probAB)
此计算得出结果 0.7530913
。
我有关节p.d.f。 我现在正在比较条件概率的理论值和我的经验值 运行 Monte Carlo 方法。 我必须重复 10,000、100,000 和 1,000,000 次抽奖。请问如何在R代码中复制?
此外,对于最后一步,条件概率, 我正在使用 Monte Carlo 方法来做到这一点。 是否有任何 R 代码可用于多元制服 分布计算条件概率?
如有任何建议,我们将不胜感激!谢谢!!
我的代码如下:
# f(x.y) = (1/4)*xy, 0<x<2, 0<y<2
# Find P(A) = P(X>1)
f <- function(x){(1/2)*x} # Marginal P(X)
probE <-integrate(f, lower = 1, upper = 2)
cat('\n Pr[ 1 < X ] is \n')
print(probE)
n <- 10000
x<-runif(n, 1,2)
probE.MC <- ((2-1)/n)*sum((1/2)*x)
cat('\n Monte Carlo Pr[1< X ] =',probE.MC,'\n')
# Find P(B) = P(Y<1)
f <- function(y){(1/2)*y} # Marginal P(Y)
probB <-integrate(f, lower = 0, upper = 1)
cat('\n Pr[ 1< Y ] is \n')
probB
typeof(probB)
n <- 10000
y<-runif(n, 0,1)
probB.MC <- ((1-0)/n)*sum((1/2)*y)
cat('\n Monte Carlo Pr[Y < 1] =',probB.MC,'\n')
# Pr[A intersect B]
# P[X>1 and Y <1]
f <- function(x,y){return((1/4)*x*y)}
n <- 100000
a11<-1; a12 <-2; a21 <- 0; a22 <-1
x <-runif(n, a11, a12)
y <- runif(n,a21, a22)
probMC <- ((a12-a11)*(a22-a21)/n)*sum(f(x,y))
probMC
typeof(probMC)
# P[A|B] = p[A intersect B]/ P(B)
probAB <- probMC/probB
首先,我重新格式化了函数并给它们单独命名。
fX <- function(x) {
0.5 * x # Marginal P(X)
}
# Find P(B) = P(Y<1)
fY <- function(y) {
0.5 * y # Marginal P(Y)
}
fXY <- function(x, y) {
1 / length(x) * sum(0.25 * x * y) # joint X,Y
}
进行多次运行的最简单方法是将 MC 代码包装在一个 for 循环中,并将每个计算保存在一个数组中。然后,最后取存储值的平均值。
所以对于 P[A] 你有:
n <- 10000
probA.MC <- numeric(n) # create the array
for (i in 1:10000) {
x<-runif(n, 1,2)
probA.MC[i] <- ((2-1) / n) * sum(0.5 * x)
}
cat('\n Monte Carlo Pr[1 < X] =',mean(probA.MC),'\n')
(我假设 probE.MC 应该是 probA.MC。)结果是 Monte Carlo Pr[1 < X] = 0.7500088
。该代码类似于 P[B],结果为 Monte Carlo Pr[Y < 1] = 0.2499819
.
我们使用 fXY 作为联合概率。
n <- 10000
probMC <- numeric(n)
for (i in 1:10000) {
x <- runif(n, 1, 2)
y <- runif(n, 0, 1)
probMC[i] <- ((a12-a11) * (a22-a21)) * fXY(x, y)
}
cat('\n Monte Carlo Pr[X,Y] =',mean(probMC),'\n')
这个结果是Monte Carlo Pr[X,Y] = 0.1882728
。
您最后一次计算应该如下所示(注意积分结果中的 probB$ 值):
# P[A|B] = p[A intersect B]/ P(B)
probAB <- mean(probMC) / probB$value
print(probAB)
此计算得出结果 0.7530913
。