Рассчитать условную вероятность совместного pdf - PullRequest
Новая
       6

Рассчитать условную вероятность совместного pdf

0 голосов
/ 21 сентября 2018

У меня есть совместный pdf И сейчас я сравниваю теоретическое значение условной вероятности и эмпирическое значение, которое я использовал для подхода Монте-Карло.Я должен сделать это в тиражах 10 000, 100 000 и 1 000 000 тиражей.Могу ли я спросить, как поместить репликации в коде R?

Также, для последнего шага, условной вероятности, я использовал подход Монте-Карло, чтобы сделать это.Есть ли какой-нибудь 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

1 Ответ

0 голосов
/ 22 сентября 2018

Сначала я переформатировал функции и дал им отдельные имена. 

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.

...