How to estimate the probability that a single value follows the same distribution as a a set of values

Question

This question is a possible duplicate of this one but I would like to go a bit further.

I have a set of values $X=x_1, x_2, \cdots, x_n$ that are iid estimates of a reference value $y$ given by a Monte-Carlo like algorithm. I know my reference value $y$ and I would like to assess the quality of the prediction given by the algorithm. An indicator of the quality of the prediction could be the answer to the following question:

What is the probability of picking $y$ using the distribution of $X$?

I know about the solution based on the rank of the reference value proposed here but it seems limited to me. In the example below the estimated values (blue bars) follow a long tailed distribution and the orange, green and red bars are three different cases of reference value. The orange and the green bars have the same rank but it seems pretty obvious to me that the probability of the orange bar being picked from the blue bars distribution is higher that for the green bar. The green and the red bars have a different rank but I would intuitively say that their probability of being picked using the blue bars distribution is quite similar.

I intuitively think there could be some kind of moment-like indicator to estimate how the distribution fits with the reference. For example:

$I=\sqrt{\frac{\sum_{i}^{n}{(x_i-y)^2}}{n}}$

But I don't know how to relate it to a probability of the reference being picked using the same distribution as the sample points.

Is there any method giving a better answer to this question than the rank-based method?

Edit: I precise that I have no assumption about the distribution of $X$.

Answer 1

There are two approaches here. For clarity, let's call your set of sampled values $X$, and your single reference value $y$.

Parametric Approach

Are the values $X$ distributed according to any known parametric distribution, for example the log-normal distribution? If so, you can estimate the parameters of this distribution ($\mu$ and $\sigma$ for the log-normal), and then use the cumulative distribution function for that distribution, with those parameters, find the p-value: the probability of finding a value as extreme or more extreme than $y$ under this distribution.

mu = 2
sigma = .5
X = rlnorm(1000, mu, sigma)
y = 25

hist(X, breaks = 20)
abline(v=y, col='red', lwd=2)

estimated.mu = mean(log(X)) # ≈ 2
estimated.sigma = sd(log(X)) # ≈ .5

# Probability of getting a value of y or lower from this distribution
p.low = plnorm(y, estimated.mu, estimated.sigma, lower.tail = T) 
# Probability of getting a value of y or higher
p.high = plnorm(y, estimated.mu, estiated.sigma, lower.tail = F)
(p.value = 2*min(p.low, p.high))
# ≈ 0.015

Non-parametric approach

If $X$ doesn't match any known distribution, you can do something a little more naive, and just count the proportion of $X$ that is above or below $y$.

# Observed proportion of values lower than y
approx.p.low = mean(y > X)
# Observed proportion of higher
approx.p.high = mean(y < X)
(p.value.alt = 2*min(approx.p.low, approx.p.high))
# ≈ 0.015

Note that neither of these approaches works properly if $X$ is bimodally distributed, since they won't pick up on cases where $y$ falls between the two modes.