5
$\begingroup$

We do a measurement on 1000 samples to detect if a chemical element A is present, and for each measurement, two cases can happen :

  • the element A is not present, and the values we get are a "background noise", with a gaussian distribution $N(100, 15)$.

  • the element A is present, and then we get values around $120$.

Given a list of 1000 measurements [101.1, 98.2, 123.12, ..., 125.5], how can we conclude with a given confidence %, that, the element A has been found at least once, or at least k times, etc. ?

Intuitively, we can look at the histogram, and see if it deviates or not from a gaussian curve. How can we turn this into a quantitative solution?

enter image description here

Are there solutions where I subtract the density of N(100, 15) (=the black curve) from the normalized histogram of the sample (=the orange curve), and see if the difference can highlight a population:

enter image description here

?

$\endgroup$

2 Answers 2

8
$\begingroup$

The natural way to model this is as a mixture of two Normal distributions. Proportion $\pi$ of the data comes from distribution 1, which has a mean of 100, a standard deviation of 15. Proportion $1 - \pi$ comes from distribution 2, which has a mean of 120, and an unknown standard deviation (unless you're able to provide one). You're primarily interested in estimating $1 - \pi$.

There are lots of tools for estimating Gaussian mixture models, but most won't quantify uncertainty in the proportions, as far as I know. I would therefore reach for Bayesian tools like Stan here. Here's a very rough implementation of this idea, based on an example from the rstan documentation. If anyone wants to edit this answer to provide more explanation, please feel free.

library(rstan)
# Simulate data
n = 1000
n_k = 50
data = c(
  rnorm(1000 - n_k, 100, 15),
  rnorm(n_k, 120, 2)
)
hist(data, breaks=101)

enter image description here

Stan model

data {
  int N;
  vector[N] y;
}
parameters {
  ordered[2] mu;
  vector<lower = 0>[2] sigma;
  // Here, Theta[1] = pi, Theta[2] = 1-pi
  simplex[2] Theta;
}
model {
  vector[2] contributions;
  // priors
  mu[1] ~ normal(100, .1);
  sigma[1] ~ normal(15, .1);
  mu[2] ~ normal(120, .1);
  sigma[2] ~ cauchy(0, 1);
  Theta ~ dirichlet(rep_vector(2.0, 2));
  
  // likelihood
  for(i in 1:N) {
    for(k in 1:2) {
      contributions[k] = log(Theta[k]) + normal_lpdf(y[i] | mu[k], sigma[k]);
    }
    target += log_sum_exp(contributions);
  }
}
stan_data = list(
  N = length(data),
  y = data
)
stan_result = rstan::sampling(mm_stan_model, data=stan_data, iter=2000)
stan_result
# Inference for Stan model: stan-c893ebbe0d5.
# 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
# mu[1]      100.00    0.00 0.10    99.80    99.94   100.00   100.07   100.20  3301 1.00
# mu[2]      120.01    0.00 0.10   119.81   119.94   120.01   120.08   120.21  3289 1.00
# sigma[1]    15.01    0.00 0.10    14.82    14.95    15.01    15.08    15.20  4030 1.00
# sigma[2]     1.97    0.03 0.75     1.15     1.56     1.82     2.14     4.03   457 1.00
# Theta[1]     0.95    0.00 0.01     0.93     0.94     0.95     0.96     0.97   910 1.00
# Theta[2]     0.05    0.00 0.01     0.03     0.04     0.05     0.06     0.07   910 1.00
# lp__     -4145.98    0.06 1.71 -4150.22 -4146.90 -4145.62 -4144.69 -4143.77   795 1.01
# 
# Samples were drawn using NUTS(diag_e) at Mon Apr 22 13:49:35 2024.
# 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).
posterior_samples = extract(stan_result)
posterior_one_minus_pi = posterior_samples$`Theta`[,2]
hist(posterior_one_minus_pi, breaks=101) # Proportion of values that are outliers

enter image description here

Update

There's actually a simpler approach here, if you are willing to assume the means and standard deviations of the two distributions, but I don't have time to go into it.

$\endgroup$
6
  • $\begingroup$ Thank you for your answer. I am not familiar with this approach. What is the conclusion after the final histogram of posterior_one_minus_pi? "The maximum likelihood is reached when $1-\pi = 0.05$", or is there another way to write the conclusion of this approach? $\endgroup$
    – Basj
    Apr 22 at 13:02
  • $\begingroup$ Sorry, yes, this won't totally make sense if you're not familiar with Bayesian methods in general. I'd suggest Doing Bayesian Data Analysis as a great introduction to the topic. You can see from the table that Theta[2] has an estimated value of 0.5, and a 95% confidence interval of [0.3, 0.7]. $\endgroup$
    – Eoin
    Apr 22 at 14:17
  • $\begingroup$ Yes @Eoin, I don't have the background, but I'll try to read about it. What does Theta (and Theta[2]) represent, and why a Dirichlet distribution? $\endgroup$
    – Basj
    Apr 22 at 14:23
  • $\begingroup$ Theta a pair of parameters indicating what proportion of the data comes from the first distribution (Theta[1]) and what proportion comes from the second (Theta[2]). $\endgroup$
    – Eoin
    Apr 22 at 14:27
  • $\begingroup$ The Dirchlet distribution is a common prior distribution for parameters that represent proportions (all of the values of theta have to add up to 1). I've used a very broad prior here, so this will have little impact on your parameter estimates. $\endgroup$
    – Eoin
    Apr 22 at 14:29
3
$\begingroup$

Answering questions such as element A has been found at least once, or at least k times, etc. is not simple. It's easier to answer questions pertaining to the whole sample rather than individual units.

What you could do is compare your obtained distribution with a reference distribution having mu and sigma equal to 100 and 15 (assuming you know these values and they are true) by using the Kolmogorov-Smirnov test. You could also use a test for checking whether your distribution is "normal" (which it wouldn't be if element A is present and you would have a bimodal distribution) by using the Shapiro-Wilk test.

You could also simply compare the mean to 100 using a (one sided?) t-test.

$\endgroup$

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.