¿Cómo usted calcular el P(x>y) para una función de densidad conjunta en R?

Question

Estoy tratando de comprender el Bayesiano AB proceso de pruebas más a fondo. Si tengo dos pruebas de forma tal que las posteriores son: $$x\sim Beta(\alpha_1, \beta_1)$$ $$y\sim Beta(\alpha_2, \beta_2)$$

Where $x$ and $s$ are posterior distributions. Then the joint pdf is of the form:

$f(x,y) = \frac{\gamma(\alpha_1+\beta_1)}{\gamma(\alpha_1)\gamma(\beta_1)}x^{\alpha_1-1}(1-x)^{\beta_1-1}\frac{\gamma(\alpha_2+\beta_2)}{\gamma(\alpha_2)\gamma(\beta_2)}y^{\alpha_2-1}(1-y)^{\beta_2-1}\quad 0\leq x \leq 1,\quad 0\leq y \leq 1$

I want to know $P(x>y)$, how do u solve the equation

$$ P(x>y) = \int_{0}^{1}\int_{y}^{1}\frac{\gamma(\alpha_1+\beta_1)}{\gamma(\alpha_1)\gamma(\beta_1)}x^{\alpha_1-1}(1-x)^{\beta_1-1}\frac{\gamma(\alpha_2+\beta_2)}{\gamma(\alpha_2)\gamma(\beta_2)}y^{\alpha_2-1}(1-y)^{\beta_2-1}dxdy$$

in R? I have the joint function coded like so

joint_dist <- function(x, y, aX, bX, aY, bY) {
  (gamma(aX + bX)/gamma(aX)*gamma(bX)*x^(aX - 1)*(1-x)^(bX-1))*(gamma(aY + 
   bY)/gamma(aY)*gamma(bY)*y^(aY-1)*(1-y)^(bY-1))
  }

How do I use the integrate function to integrate from y to 1 and then 0 to 1? Also I'm aware of Evan Miller's summation derivation, but I'm under the impression that he uses a $Beta(1,1)$ antes para el control y la variación. Me gustaría ser capaz de elegir una diferente previo para cada uno. Tal vez un uniformative antes de la variación, pero más informado previo para el control. Tal vez estoy equivocado aquí.

Answer 1

Cómo sobre el uso de MonteCarlo a la integración?

#You can put your value here. 
aX=3 ;bX=7; aY=3; bY=7

#MC Integration
n_sim=1000000   # number of simulation
u<-runif(n_sim) # 
C<-(1-pbeta(u,aX,bX))*dbeta(u,aY,bY)
m<-mean(C)  #integration result
s<-sd(C)    # standard deviation of the result.

Usted puede averiguar el intervalo de confianza del cálculo utilizando la desviación estándar. Por ejemplo, el IC del 95% sería (m-s/sqrt(n_sim)*1.96,m+s*1.96/sqrt(n_sim)).

Cuando me puse a X y y tienen la misma distribución, tengo el 50% como resultado.

Podría ser la mejor manera de calcular y tal vez alguien podría publicar como otra respuesta.

Yo no soy un hablante nativo de inglés por lo que cualquier revisión de esta respuesta se agradece. Gracias de antemano!

edit1 Me gustaría añadir algo más de detalle.

Digamos

$1-F_{\alpha_1\beta_1}(y)=\int_{y}^{1}\frac{\gamma(\alpha_1+\beta_1)}{\gamma(\alpha_1)\gamma(\beta_1)}x^{\alpha_1-1}(1-x)^{\beta_1-1}dx$

y

$f_{\alpha_2\beta_2}(y)=\frac{\gamma(\alpha_2+\beta_2)}{\gamma(\alpha_2)\gamma(\beta_2)}y^{\alpha_2-1}(1-y)^{\beta_2-1}$

a continuación, $P(x>y) =\int_{0}^{1} (1-F_{\alpha_1\beta_1}(y))f_{\alpha_2\beta_2}(y) dy$.

Y esta expresión es en realidad llegar expectativa $E[(1-F_{\alpha_1\beta_1}(y))f_{\alpha_2\beta_2}(y)]$ cuando y tiene una distribución uniforme.

Así que lo que hice es esencialmente la simulación de esta expectativa. Tomé 1000000 (n_sim) muestras de la distribución uniforme y obtenga $(1-F_{\alpha_1\beta_1}(y))f_{\alpha_2\beta_2}(y)$ por cada cada una de las muestras (que es C) y calcula la media y la desviación estándar. Esto significa que es imparcial estimador de la integración. Gracias al teorema del límite central, la media de C tiene una distribución asintótica normal y esta es la razón por la que podemos obtener el intervalo de confianza también.

Además, supongo que se puede interpretar la integración como $E[(1-F_{\alpha_1\beta_1}(y))]$ cuando y tiene una distribución beta con $\alpha_2$, $\beta_2$.

y<-rbeta(n_sim,aY,bY)
C2<-(1-pbeta(y,aX,bX))
mean(C2)

Esto daría a la misma media. Pero en este caso la desviación estándar parece diferente. No he averiguado por qué es lo que tengo que hacer más investigación sobre esta diferencia.

Sin embargo, el método de obtener el intervalo de confianza usando el uniforme de la distribución de referencia.