Fit probability distribution from probabilistic observations?

Question

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:

https://pastebin.com/xcUdS7ch

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

Answer 1

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 *)
SeedRandom[12345]
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.

Answer 2

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.