# Detect rare high-value measurements in a series of measurements

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?

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:

?

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)


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  ## 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. • 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?
– Basj
Apr 22 at 13:02
• 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].
– Eoin
Apr 22 at 14:17
• 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?
– Basj
Apr 22 at 14:23
• 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]).
– Eoin
Apr 22 at 14:27
• 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.
– Eoin
Apr 22 at 14:29

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.