Trying to calculate the probability that one distribution (fitted to data) is larger than another one. I've tried this:

Probability[x < y, {x \[Distributed] PoissonDistribution[26.],
   y \[Distributed] SkewNormalDistribution[11.7, 7.93, 3.32]}]

This just crashes kernel. Also tried:

NProbability[x < y, {x \[Distributed] PoissonDistribution[26.],
  y \[Distributed]    SkewNormalDistribution[11.7, 7.93, 3.32]}]

This just gives:

NProbability[x < y, {x, y} \[Distributed] ProductDistribution[PoissonDistribution[26.],
 SkewNormalDistribution[11.6727, 7.93249, 3.31684]]]

What am I doing wrong?

Thanks very much!

  • $\begingroup$ Your comparison is not well defined for one variable from a discrete distribution and another from a continuous one. $\endgroup$ Jun 9, 2017 at 18:11
  • 1
    $\begingroup$ @DavidG.Stork I think the comparison is well defined although certainly not usual. $\endgroup$
    – JimB
    Jun 9, 2017 at 18:19
  • $\begingroup$ @JimBaldwin, any further thoughts on why Probability does not give an answer even though the problem is well-defined? $\endgroup$
    – rfrasier
    Jun 12, 2017 at 1:52
  • $\begingroup$ @rfrasier More thought and trying things with Probability but little insight. I tried what I thought would be a simpler example using what was stated in the "Details and Options" of Probability but the result did not make things clear: Expectation[ Probability[x < y, y \[Distributed] BinomialDistribution[26, 8/10], Assumptions -> 0 <= x <= 20 && x \[Element] Integers], x \[Distributed] Binomial[20, 9/10]]. I say "well defined" in a mathematical sense but that doesn't make it work in Mathematica. $\endgroup$
    – JimB
    Jun 12, 2017 at 3:45
  • $\begingroup$ @JimBaldwin, thanks for your follow-up. I took your meaning of "well defined" to be the mathematical one. :-) Just noticed a small typo in your code---I believe you meant x is distributed according to BinomialDistribution[20, 9/10], in which case, your code evaluates to a Rational number with a Real value of approx. (*0.829314*). All of our solutions are still numerical approximations without any form of the series form. I would like to investigate further to see if there is an analytic form. $\endgroup$
    – rfrasier
    Jun 15, 2017 at 4:06

3 Answers 3


Here's a brute-force approach by sampling:

nsim = 10000000;
x = RandomVariate[PoissonDistribution[26], nsim];
y = RandomVariate[SkewNormalDistribution[11.7, 7.93, 3.32], nsim];
N[Count[x - y, u_ /; u < 0]/nsim]
(* 0.1270612 *)

(This, of course, assumes the random variables are independent and it does not at all account for the fitting of the skew normal distribution.)

A more direct approach is conditioning on the values of the Poisson distribution:

Sum[PDF[PoissonDistribution[26], x] (1 - CDF[SkewNormalDistribution[11.7, 7.93, 3.32], x]),
  {x, 0, 40}]
(* 0.127095 *)

It looks like no more than 40 terms are needed.

Now if the error in the estimation of the parameters of the skew normal distribution needs to be taken into account (i.e., getting reasonable confidence intervals for the probability that X < Y), then a bootstrap process could be considered either by re-sampling the data or sampling from the estimated skew normal distribution.


If the underlying problem is to estimate $Pr(X<Y)$ from a known Poisson distribution with mean 26 and from a skew normal distribution with parameters estimated by a sample, then all of the methods mentioned in all of the answers work fine. However, the golden rule of estimation is to always obtain a measure of precision along with the estimate. (That's well known in most fields but not always practiced.)

One way to obtain an estimate of precision is to use some sort of bootstrap process. Below is one such approach. (It is written in a step-by-step approach which is not necessarily computationally efficient. But that can come later.)

(* Suppose we have a random sample of size n from a skew normal distribution *)
n = 500;
x = RandomVariate[SkewNormalDistribution[12, 8, 3], n];

(* Find maximum likelihood estimates *)
mle = FindDistributionParameters[x, SkewNormalDistribution[μ, σ, α]]
(* {μ -> 12.311751931399243, σ -> 7.436365977762536, α -> 2.588196061903045} *)

(* Get parametric bootstrap samples of size 100 from the estimated skew normal distribution *)
nboot = 1000;
bootMLE = Table[{μb, σb, αb} /. 
  FindDistributionParameters[RandomVariate[SkewNormalDistribution[μ, σ, α] /. mle, 100], 
    SkewNormalDistribution[μb, σb, αb], {{μb, μ /. mle}, {σb, σ /. mle}, {αb, α /. 
       mle}}], {i, nboot}];

(* Determine Pr(X<Y) for each bootstrap sample *)
probBoot = 
  Table[Sum[PDF[PoissonDistribution[26], x] 
    (1 - CDF[SkewNormalDistribution[bootMLE[[i, 1]], bootMLE[[i, 2]], 
         bootMLE[[i, 3]]], x]), {x, 0, 40}], {i, nboot}];

(* We can summarize our knowledge about the estimate of the probability in a variety of ways *)
(* Standard error *)
probBootSE = StandardDeviation[probBoot]
(* 0.0216423 *)

(* Approximate 95% confidence intervals *)
prob95ConfInt = 
 Sort[probBoot][[{Floor[0.025 nboot], Floor[0.975 nboot] + 1}]]
(* {0.113747, 0.197672} *)

(1,000 bootstraps takes about 2.5 minutes to run.)

  • $\begingroup$ answered before reading yours, happy to delete mine if you desire - same general idea as your final... $\endgroup$
    – ciao
    Jun 10, 2017 at 1:15
  • $\begingroup$ @ciao Yours is good. No reason to delete. $\endgroup$
    – JimB
    Jun 10, 2017 at 2:33

This will get you your result:

Tr[PDF[PoissonDistribution[26.], #]*
   SurvivalFunction[SkewNormalDistribution[11.7, 7.93, 3.32], #] &@Range[0, 40]]


Range[0,40] suffices for 6 significant figures, compare with exact parameters:

Tr[PDF[PoissonDistribution[26], #]*
     SurvivalFunction[SkewNormalDistribution[117/10, 793/100, 332/100], #] &@ Range[0, 40]] // N[#, 6] &



I was able to get the solution to work when I changed the Method for the NProbability function to "MonteCarlo". Here's the code I used to arrive at the answer:

NProbability[x < y, 
 {x \[Distributed] PoissonDistribution[26.], 
  y \[Distributed] SkewNormalDistribution[11.7, 7.93, 3.32]}, 
 Method -> {"MonteCarlo", PrecisionGoal -> 6}]



I am unsure why using the Probability function to get the answer. Perhaps someone else could address this related issue.

  • $\begingroup$ The documentation does not require && but only indicates that you can use && instead of a List. Change that example in the documentation to a list and it works fine. (But the wording in the documentation could make that clearer.) Also while one could discretize the continuous distribution, it is not a requirement. $\endgroup$
    – JimB
    Jun 10, 2017 at 14:14
  • $\begingroup$ @JimBaldwin, noted with thanks! I wil edit my answer to reflect. $\endgroup$
    – rfrasier
    Jun 12, 2017 at 1:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.