2
$\begingroup$

I am trying to solve the following problem, but I cannot find google-relevant keywords for it (I'm not a statistician) :

  • In a virtually infinite population
  • All individues can have trait $A$ with a known probability $p_A$
  • Individues who have no trait $A$, have an unknown probability $p_B$ of having trait $B$
  • Individues can have either trait A, or trait B, or none. They cannot have two traits.
  • In a sample of size $N$, we observe experimentally $X = #A + #B$, where $#A$ is the number of individues in the sample that have trait $A$ and $#B$ is the number of individues in the sample that have trait $B$. #A and #B are latent variables, we cannot observe them directly. Only their sum can be observed.

Given that The prior on $p_B$ is uniform in [0,1], $p_A$ is known. and X is known (but not #A nor #B separately) : What is the updated posterior probability distribution of $p_B$ ?

I am looking for a solution that would still be correct even with $p_A \cdot N \leq 5$ and $p_B \cdot N \leq 5$

I have the feeling that this should have already been solved in the literature, but I can't find it, nor the relevant keywords to find it

Intuitive extreme example:

  1. I have 1000 individues
  2. I have $p_A = 10^{-7}$
  3. I observe experimentally X = 500
  4. From this observation of X among N, I can deduce with an important confidence that $p_B \approx 0.5$
$\endgroup$
5
  • 1
    $\begingroup$ Can people have both traits? $\endgroup$
    – bdeonovic
    Oct 29, 2021 at 2:14
  • 1
    $\begingroup$ what do you mean by the sum of $A$ and $B$? do you mean in a sample $N$ you know how many have either trait $A$ or $B$? $\endgroup$
    – bdeonovic
    Oct 29, 2021 at 2:18
  • $\begingroup$ in your sample do you know how many are trait $A$, or only the total with either trait? $\endgroup$
    – bdeonovic
    Oct 29, 2021 at 2:21
  • $\begingroup$ Hi @bdeonovic thanks for your interest. I updated the post to be more precise. People cannot have two traits. I only know the sample size $N$, the total number $X$ of people with either traits in the sample, and $p_A$. And I suppose a uniform(0,1) prior for $p_B$ $\endgroup$ Oct 29, 2021 at 3:57
  • $\begingroup$ You ask for keywords. In the literature on contingency tables there is the problem of misclassification where you know which row the observation is in but not the column. That is not exactly your situation but perhaps might give you some leads. Try putting contingency table misclassification into your favourite search engine. $\endgroup$
    – mdewey
    Oct 29, 2021 at 13:12

1 Answer 1

7
$\begingroup$

I have no idea what such a problem corresponds to in the literature, but the problem as you have stated it is quite simple to do inference on:

You have

$$ \begin{align*} X|p_b &\sim \text{Binomial}(N,p_a + (1-p_a)p_b)\\ p_b &\sim \text{Beta}(1,1) \end{align*} $$

The posterior is then proportional to $$ \begin{align*} p(p_b|X) &\propto p(X|p_b)p(p_b)\\ &\propto (p_a + (1-p_a)p_b)^X(1 - (p_a + (1-p_a)p_b))^{N-X} \end{align*} $$

This doesn't correspond or simplify to any named distribution (as far as I know). However we can easily sample from the posterior to do inference on $p_b$. For example in R you can do this using rstan library

library(rstan)

model<-"
data {
  int<lower=0> X;          
  int<lower=0> N;          
  real<lower=0,upper=1> pA;
}
parameters {
  real<lower=0,upper=1> pB;
}
model {
  target += binomial_lpmf(X | N, pA + (1-pA) * pB);
  target += beta_lpdf(pB | 1.0, 1.0);
}
"

Lets test this for $N=1000$, $X=500$, and $p_a=0.25$

data <- list(
  pA=0.25,
  X=500,
  N=1000
)

fit <- stan(model_code=model, data=data)

We get

Inference for Stan model: 54c52a299004e265657393fabf67c2e0.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

      mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
pB    0.33    0.00 0.02  0.29  0.32  0.33  0.35  0.38  1209    1
lp__ -5.68    0.02 0.71 -7.76 -5.81 -5.40 -5.23 -5.18  1498    1

Samples were drawn using NUTS(diag_e) at Fri Oct 29 07:48:30 2021.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1)

so the posterior mean is $0.33$ with a 95% credible interval $(0.29,0.38)$.

EDIT: As @whuber points out if $U=(1-p_A)(1-p_B)$ then $$U|X \sim \text{Beta}(N-X+1,X+1)$$ The inverse transform gives $p_B = 1 - \dfrac{U}{1-p_A}$. As @whuber points out you can do many things now just using standard software for beta distributions like sample from the beta and then apply the inverse transform to get samples from posterior of $p_B$. Since the inverse transform is affine you can easily compute expectations $$ \text{E}[p_B|X] = 1 - \dfrac{\text{E}[U|X]}{1-p_A} = 1 - \dfrac{\tfrac{N-X+1}{N-X+1 + X+1}}{1-p_A} $$ which is equal to $1/3$ for the example I used.

$\endgroup$
6
  • $\begingroup$ Let $U = (1-p_b)(1-p_a).$ Then, conditional on $X,$ $U$ has a Beta$(N-X+1, X+1)$ distribution. $\endgroup$
    – whuber
    Oct 29, 2021 at 18:16
  • $\begingroup$ @whuber Do you think it wouldt be possible to use this to find a faster way to compute the posterior ? $\endgroup$ Oct 29, 2021 at 19:13
  • $\begingroup$ It depends on what might mean by "compute." But knowing $p_a=1-U/(1-p_b)$ is related so simply to a Beta distribution means you can easily and efficiently calculate most of its properties as well as draw random values from it using standard functions supplied by most statistical software. $\endgroup$
    – whuber
    Oct 29, 2021 at 19:22
  • 1
    $\begingroup$ @Sextus I don't think so. Notice the parameters are $N-X+1$ and $X+1,$ in that order. $\endgroup$
    – whuber
    Oct 29, 2021 at 21:17
  • 1
    $\begingroup$ Ah I see if $U$ is beta distributed then so is $1-U$ $\endgroup$ Oct 29, 2021 at 21:19

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.