5
$\begingroup$

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!

| improve this question | |
$\endgroup$
  • $\begingroup$ Your comparison is not well defined for one variable from a discrete distribution and another from a continuous one. $\endgroup$ – David G. Stork Jun 9 '17 at 18:11
  • 1
    $\begingroup$ @DavidG.Stork I think the comparison is well defined although certainly not usual. $\endgroup$ – JimB Jun 9 '17 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 '17 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 '17 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 '17 at 4:06
6
$\begingroup$

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.

Update

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.)

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

This will get you your result:

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

0.127095

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] &

0.127095

| improve this answer | |
$\endgroup$
2
$\begingroup$

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}]

returns

(*0.126862*)

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

| improve this answer | |
$\endgroup$
  • $\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 '17 at 14:14
  • $\begingroup$ @JimBaldwin, noted with thanks! I wil edit my answer to reflect. $\endgroup$ – rfrasier Jun 12 '17 at 1:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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