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.)
Probability
does not give an answer even though the problem is well-defined? $\endgroup$ – rfrasier Jun 12 '17 at 1:52Probability
but little insight. I tried what I thought would be a simpler example using what was stated in the "Details and Options" ofProbability
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:45x
is distributed according toBinomialDistribution[20, 9/10]
, in which case, your code evaluates to aRational
number with aReal
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