1

I would like to compute the probability $P(X > Y)$ with R. I used JAGS to sample from the posterior distribution of each variable, so I have a Markov chain for each variable (of length $3\times 10^{5}$).

Could someone point me in the right direction?

EDIT: preferably with a Bayesian mindset.

FatherNucleus
  • 21
  • 4
  • 3
    Not worth a whole answer but just count the number of pairs where $X>Y$ and divide by the total number of pairs. More formally, you are using the MCMC chain to approximate the expected value of $\mathbb{I}(X>Y)$. – deasmhumnha Jun 12 '18 at 22:03
  • @Dezmond: That is absolutely worth a whole answer. It is the only possible thing that can be done on the basis of the information in the question. – Ben Jun 13 '18 at 00:11
  • @FatherNucleus: You say you have simulated the posterior, but the probability you have specified is not expressed to be conditional on any past observations. If you have correctly specified your question then what you are asking for is impossible, since you are asking for a prior predictive probability, when your only information is a simulated posterior sample. Please check the framing of your question - did you intend for the probability of interest to be a posterior prediction? If so, please adjust your question accordingly. – Ben Jun 13 '18 at 00:19
  • @Ben: Perhaps my question was a little unclear. I'll give a little more detail here, but I think the question contains all the information that can be worked with. X and Y are two parameters in a logit regression. I've used uninformative priors for them and obtained two Markov chains accordingly through JAGS. I agree with you that Dezmond's answer seems like the only thing to do. – FatherNucleus Jun 13 '18 at 02:24
  • @DezmondGoff: I would like to accept your comment as answer. This seems like the only sensible thing to do, and in fact, was the sort of thing I was looking for. – FatherNucleus Jun 13 '18 at 02:25
  • I think a mod would have to migrate it, right? – deasmhumnha Jun 13 '18 at 04:12

1 Answers1

1

Based on your question and your elaboration in the comment section, it appears that what you are trying to find is the posterior probability that one parameter $X$ exceeds another parameter $Y$. (It is unusual to represent parameters with English letters; usually we represent data with English letters and parameters with Greek letters.)

You report have a posterior simulation using MCMC methods in JAGS, and you sampled from the posterior distribution of each variable. Be careful here: parameters are usually not independent a posteriori, so you need to sample from a chain that generates the joint posterior of the parameters; don't sample separately from two chains generating the marginal posteriors. Assuming you have done this correctly, the simplest way to estimate the probability is via the proportion of cases where the posterior event of interest holds in your simulation. (Hat-tip to Dezmond Goff in the comments for being the first to point this out.)

Suppose you have $b$ burn-in iterations in your chain, and then an additional $n$ simulated values, so your simulated values are $\{x_k, y_k\}_{k=1}^{b+n}$. Using these simulated MCMC values you can estimate the posterior probability of interest by discarding the burn-in simulations and using the rest:

$$\hat{\mathbb{P}}(X>Y|\text{Data}) = \frac{1}{n} \sum_{k=b+1}^{b+n} \mathbb{I}(x_k > y_k).$$

Without any additional information about the distributions, this is the best you can do. You are using your simulated MCMC values to directly estimate the posterior probability of your event. Ergodic properties of MCMC methods and the laws-of-large numbers ensure that this estimator converges to the rue probability as $n \rightarrow \infty$, so you should get a good answer with $n=3 \times 10^5$.

Ben
  • 91,027
  • 3
  • 150
  • 376
  • Thank you for your answer. Just out of curiosity, how would something like this do: $\mathrm{P}(X > Y|\mathrm{Data}) \approx \frac{1}{n^2} \sum_{j=1}^{n} \sum_{i=1}^{n} \mathrm{I}(x_i > y_j)$? – FatherNucleus Jun 13 '18 at 12:47
  • That is no good - you would be matching $x$ and $y$ values from different iterations of the simulation, which is equivalent to treating them as independent. – Ben Jun 13 '18 at 12:50