From unevenly spaced values, a normal distribution that sums to the mean

Question

I have input data that is $N$ items and each has several properties, one of which is called, say, $P$. Their values of $P$ are unevenly spaced between 0 and 1, and may even repeat between items. For example, suppose $N=8$ and the list of values for $P$ looks like this:

$$[0.1, 0.2, 0.2, 0.3, 0.4, 0.5, 0.8, 0.9]$$

(Recall, even though 0.2 appears twice in that list, they correspond to different items that differ in other properties).

From these items, I want to draw $M$ samples, and $M >> N$ so naturally there will be repeated draws of the same item. Let's say $M=100$.

But, crucially, I want the distribution of $P$ among the drawn samples $M$ to be approximately normal, having a specific target mean $\mu_P$ and standard deviation $\sigma_P$ value. To again pick a specific example, let's say I want $\mu_P = 0.33$ and $\sigma_P = 0.2$. Keep in mind that these are my targeted values of $\mu$ and $\sigma$, they are NOT the natural mean and stdev of the original list of values $P$.

Here's what I've figured out so far: In the simpler case where the values of $P$ are spaced out evenly between 0 and 1 with no repeat values, this would be relatively straightforward. I would compute the Gaussian equation with $\mu = 0.33$ and $\sigma = 0.2$ for each value of $P$. Then, since $P$ spans a finite range and $\mu$ is not in the dead center of it, I would also need to cut off any $P$ values above $\mu+(\mu-P_{min})$, or 0.66 in this case - otherwise, the higher the $\sigma$, the higher the error that would be introduced by drawing from the extreme right-tail values when there are no correspondingly extreme left-tail values to draw from. (If $\mu$ were instead in the upper half of the range of $P$, I would do the same cutoff on the left tail instead). Then I would normalize the resulting values so that they sum to 100, and draw that many times from each. In other words, the number of times $m$ to draw item $i$ would be:

$$g_i = (e^{-(P_i-\mu_P)^2/2\sigma_P^2} )/ \sqrt{2\pi\sigma_P^2}$$ $$g_i(P_i>0.66) = 0$$ $$m_i = 100\frac{g_i}{\sum{g_i}}$$

Obviously we will have to round these to integers and so we won't hit our target $\mu$ exactly, but the higher the number of draws $M$ gets, the closer we'll get. That's fine. Again, that all works (I've tested it in Python) if they were spaced out evenly in $P$ and there were no repeats.

If we add repeated $P$ values into the mix, it's not so hard: for item $i$, if its value of $P_i$ repeats twice, I just divide each one's number-of-times-to-draw $m_i$ by 2 so that that value of $P$ doesn't get oversampled.

But what if it's not quite an exact repeat but rather two very very close values, or in general, what if $P$ is not evenly spaced? Or what if there is both uneven spacing and repeat values, as in the very first example list of $P$ values I gave? I'm trying to come up with a more universal solution, and not sure where to go. I am doing this in Python, so if a solution isn't particularly elegant or analytical, that's okay (i.e. numerical/programming approaches are okay). Any tips would be appreciated! Feel free to point me to further reading, I searched quite a bit to no avail but I may be using the wrong search terms.

Edit: This is different from flagged possible duplicate "How to simulate data that satisfy specific constraints such as having specific mean and standard deviation?." While I found that question a very interesting read, unless I'm misunderstanding it, it seems that that question is focused on generating new, simulated data that matches a certain constraint of $\mu$ and $\sigma$. In my case here, I am stuck with my $N$ items and their corresponding $P$ values. They are real items and I can't make new ones up, nor can I change the value of their property $P$. From these items, I need to determine how many times to pick each item $i$ such that my final dataset, in aggregate, has (close to) the desired $\mu$ and $\sigma$.

Edit2: Commenter @whuber has asked for more information on the problem I'm trying to solve. The problem I am trying to solve is this: I have a set of $N$ blocks of simulated material (which in my post I generically called "items"). Each block has various properties that I can measure - thermal conductivity, electrical conductivity, porosity, etc. I specified bounds of 0 to 1 for simplicity in this example, but the span could be different - either way there is a finite span. I want to choose some number $M>>N$ of these blocks (and I'm allowed to pick the same one more than once) such that taken in aggregate, their porosity averages out to appx. some value $\mu$, but among the blocks' various values of porosity there is approximately a standard deviation of $\sigma$.

Answer 1

Analysis

Let's begin with some observations.

1. You have a great deal of flexibility, because reproducing the mean and standard deviation imposes only two constraints. A third constraint is that the sum of the selection probabilities must equal unity. Thus, with more than three data values you should expect there to be an entire manifold of solutions.

2. "Repeats" of values are meaningless because the question concerns the distribution of the random variable defined by these values. Thus, there will be nothing we can do to specify how often to select the two items with common values of $0.20$. You may select them with arbitrary probabilities that sum to the required probability associated with $0.20$.

The request for a Normal-like distribution is interesting because of one characterization of Normal distributions: the maximum-entropy distribution with a given mean and standard deviation and supported on a given interval is always a truncated Normal distribution. Let's use this as a point of departure in formulating the problem.

Given $n \ge 3$ distinct values $x_1, x_2, \ldots, x_n$, a mean $\mu$, and a standard deviation $\sigma$, find (if possible) a maximum-entropy probability distribution $\mathbf{p} = p_1, p_2, \ldots, p_n$ (meaning the chance of $x_i$ is $p_i$) for which the mean of $\mathbf{p}$ is $\mu$ and its standard deviation is $\sigma$.

This optimization problem translates to the following:

$$\text{Maximize } \sum_{i=1}^n -p_i \log(p_i)$$

$\text{subject to}$

$$\eqalign{ &p_i\ge 0, i=1, 2, \ldots, n; \\ &\sum_{i=1}^n p_i=1; \\ &\sum_{i=1}^n p_i x_i = \mu; \\ &\sum_{i=1}^n p_i x_i^2 = \mu^2 + \sigma^2. }$$

Solution

Introducing Lagrange multipliers produces the solution

$$p_i = C \exp\left(\lambda_1 x_i + \lambda_2 x_i^2\right)$$

where the $\lambda_j$ are the Lagrange multipliers associated with the last two constraints. $C$ is a normalization constant (enforcing the sum-to-unity constraint). This implies the Lagrange multipliers satisfy the system of (nonlinear) equations

$$\eqalign{ 0 = & \sum_{i=1}^n (x_i - \mu) \exp(\lambda_1 x_i + \lambda_2 x_i^2) \\ 0 = & \sum_{i=1}^n (x_i^2 - (\mu^2 + \sigma^2)) \exp(\lambda_1 x_i + \lambda_2 x_i^2). }$$

There are many ways to solve these equations. One, for instance, is simply to square both left hand sides and minimize their sum: a solution is acceptable when that sum is sufficiently close to zero. To four significant figures, the solution to the example given in the question with $x=0.1, 0.2, 0.3, 0.4, 0.5, 0.8, 0.9$, $\mu=0.33$, and $\sigma=0.2$, is

$$p_1, \ldots, p_n = 0.1965, 0.2049, 0.1975, 0.1759, 0.1447, 0.05029, 0.03019.$$

Comments

It's straightforward to show there can exist a distribution with the given values of $x_i, \mu$, and $\sigma$, if and only if (a) the intended mean $\mu$ lies within the range of the $x_i$ and (b) the intended variance $\sigma^2$ is no greater than $(\max(x_i)-\mu)(\mu-\min(x_i))$. This latter is the largest possible variance of any random variable with the same range as the $x_i$ and mean $\mu$.

The next two comments relate insights afforded by the analytical form of the solution--something that would be unavailable to us if we merely were to dump the original "maximize the entropy subject to these constraints" problem into a black-box optimizer.

The formula $p_i = C \exp\left(\lambda_1 x_i + \lambda_2 x_i^2\right)$ looks just like that of the density of a Normal distribution, as intended. Do not confuse it with a Normal distribution: these are just a finite number of discrete probabilities whose graph (as in the figure) happens to fall along a Gaussian curve. (In case the $x_i$ form an arithmetic progression this distribution will indeed approximate a discretized version of a truncated Normal distribution.)

There is no guarantee that the $p_i$ will taper to lower values away from a central mode. If you want to enforce that, you will have to give up on at least one of the constraints (on the mean and standard deviation) and impose instead the constraint $\lambda_2 \le 0$. As an example, here is the solution to the problem in the question with $\sigma=0.3$ instead of $0.2$:

Requiring a relatively large standard deviation has forced most of the probability out to the extreme values of the $x_i$.

---

Code

Here is R code used to implement a solution and plot the results.

x <- c(0.1, 0.2, 0.2, 0.3, 0.4, 0.5, 0.8, 0.9) # Values of any random variable
mu <- 0.33                                     # Intended Mean
sigma <- 0.3                                   # Intended SD (positive)
tolerance <- 1e-7                              # Desired relative accuracy
#
# Check feasibility.
#
v <- -prod(range(x) - mu)
if (v < 0)
  stop("Impossible: intended mean ", mu, " is not in the range of the data.")
if (v < sigma^2)
  stop("Impossible: intended SD ", sigma, 
       " is inconsistent with the data and mean of ", mu, ".")
#
# Make sure values in `x` are unique.
#
x <- unique(x)
#
# Define the objective function.
#
f <- function(theta) {
  p <- exp(theta[1]*x + theta[2]*x^2)
  p <- p / sum(p)
  sum((x - mu)*p)^2 + sum((x^2 - (mu^2 + sigma^2))*p)^2
}
#
# Find a solution.
#
fit <- optim(c(0, 0), f, control=list(reltol=1e-15))
if (zapsmall(c(fit$value, mu^2, sigma^2))[1] != 0) 
  warning("Could not find a solution.")
theta.hat <- fit$par
p <- exp(theta.hat[1]*x + theta.hat[2]*x^2)
p <- p / sum(p)
#
# Double-check.
#
mu.hat <- sum(p*x)
sigma.hat <- sqrt(sum(p*x^2) - sum(p*x)^2)
message("Sum: ", sum(p), "; mean: ", signif(mu.hat, 4), "; 
        SD: ", signif(sigma.hat, 4), ".")
message("Entropy: ", signif(-sum(p*log(p)), 4))
if (abs(mu.hat - mu) > tolerance*(abs(mu)+abs(mu.hat)) |
    log(sigma.hat / sigma) > tolerance)
  stop("Unable to find an accurate solution.")
#
# Display the solution.
#
plot(x, p, ylim=c(0, max(p)), type="n", main="Maximum-Entropy Solution",
     sub=paste("sigma =", signif(sigma, 3)))
invisible(mapply(function(x,p) lines(c(x,x),c(0,p), col="Gray", lwd=2), x, p))
points(x, p, pch=19)