Consider an observable $X$ which follows some unknown statistical distribution $P$, with elements $x\in X$ in the interval $-1<x<1$.

Consider also an observable $Y$ which follows a uniform random distribution, with elements $y\in Y$ in the interval $-1<y<1$.

Let's say we make a sequence of observations $O$, which consist of elements $(z,p)$, such that:

  • $z\in X$ with probability $p$
  • $z\in Y$ with probability $1-p$

Is there a way to use Mathematica to obtain some sort of "best fit" for the probability distribution $P$ of observable $X$ from the "$Y$ contaminated" observations $O$?

For example, given the particular sequence of observations O:


what is a best fit for probability distribution $P$?

  • 1
    $\begingroup$ How does $Y$ have a normal distribution but is restricted to an interval between -1 and +1 ? $\endgroup$
    – JimB
    Jan 9, 2020 at 22:43
  • 1
    $\begingroup$ @JimB essentially, I mean y=2 RandomReal[]-1. Uniformly drawn random numbers in this interval. I'll change the term to uniform random distribution. $\endgroup$
    – Kagaratsch
    Jan 9, 2020 at 23:21
  • 1
    $\begingroup$ You might also initially consider relaxing the "unknown" distribution to a "known" distribution (at least up to the specific parameters) and solve that problem first. $\endgroup$
    – JimB
    Jan 9, 2020 at 23:26
  • $\begingroup$ @JimB oh, yes, that sounds interesting, as well! Unfortunately, I'm not sure where to start. I know there are distribution fitting functions in mathematica, but don't see a clear way to write syntax that would treat such a particular problem... $\endgroup$
    – Kagaratsch
    Jan 9, 2020 at 23:29

2 Answers 2


This is an extended comment. What you have described for each $z$ is a mixture distribution and with a different but known probability for each observation. If you know the form for each of the two distributions (i.e., everything known about the distribution except the specific parameters), then you could use maximum likelihood to estimate the parameters. (And you can get estimates of precision for the parameter estimates but I won't deal with that at this point.)

Below I assume that the two distributions are beta distributions. Suppose your data observations are stored in o:

(* Define the mixture distribution *)
mixture[pr_?NumericQ, a1_?NumericQ, b1_?NumericQ, a2_?NumericQ, b2_?NumericQ] :=
  MixtureDistribution[{pr, 1 - pr}, {BetaDistribution[a1, b1], BetaDistribution[a2, b2]}]

(* Get the known probabilities from your data *)
p = o[[All, 2]];

(* Generate random sample from a mixture distribution *)
z = RandomVariate[mixture[#, 2, 6, 1, 1], 1] & /@ p // Flatten;

(* Find maximum likelihood estimates *)
mle = FindMaximum[{Sum[LogLikelihood[mixture[p[[i]], a1, b1, a2, b2], {z[[i]]}], {i,  Length[z]}], 
   a1 > 0 && b1 > 0 && a2 > 0 && b2 > 0}, {{a1, 2}, {b1, 6}, {a2, 1}, {b2, 1}}]
(* {170.911, {a1 -> 1.99437, b1 -> 5.96931, a2 -> 1.00444, b2 -> 1.02215}} *)

I not sure there's much you can do if one of the distributions is "unknown". In that case there would seem to be an identifiability issue. One might get a hint of the shape of the unknown distribution by looking only at the observations where the probability was very close to 1 (although that would limit considerably the number of observations available).

You might want to consider consulting with a statistician. No need to re-invent the wheel.


Not a super great answer, but you can extract a rough distribution $P$ by removing the data points related to $Y$ as follows.

First, bin your data. In this case, put it into 20 unique bins from $-1$ to $1$.

f[{x_, p_}] := Floor@(10 x);

bins = GatherBy[Sort@dat, f];

counts = Length /@ bins;
(* {22, 22, 19, 17, 25, 27, 32, 56, 62, 222, 213, 62, 38, 31, 34, 32, 33, 24, 14, 15} *)

You can average the probabilities for each bin in order to get the probability that the points in the bin are associated with $X$

probs = Last@Transpose@(Mean /@ bins);
(* {0.467087, 0.540432, 0.466568,..., etc. }

Now multiply the counts in each bin by the probability that the point is from $X$. Note that the distribution on $Y$ is not used or needed.

xcounts = probs counts
(* {10.2759, 11.8895, 8.86479, 8.71304, 12.7415, 11.421, 16.5942, \
    25.6964, 32.5254, 106.204, 105.917, 33.3226, 13.9169, 14.308, \
    18.7334, 17.9122, 14.0052, 13.1306, 4.85573, 8.04254}

These counts of $X$ in the bins are centered on values of $X$ of

Range[20]/10 - 1.05

(* {-0.95, -0.85, -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, \
    -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95} *)

From here you can use standard methods of inferring a distribution of $X$. It could possibly help you identify the shape of the distribution, and then you could use @JimB's method to finish it off.

  • 2
    $\begingroup$ I think this answer is better than you suggest. Your approach doesn't require a known distribution form and could be modified to allow for a nonparametric density estimate based on the estimated count (xcounts) from distribution $P$. (Is such an approach in the literature? I don't know. Certainly there are papers about nonparametric estimates of mixture distributions but maybe not where the mixture proportions are known and vary for each observation.) $\endgroup$
    – JimB
    Jan 10, 2020 at 16:41
  • $\begingroup$ It's an unusual, interesting problem for sure. $\endgroup$
    – MikeY
    Jan 10, 2020 at 20:27

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.