Sampling data to have specific mean and standard deviation [duplicate]

Question

I have a data that I want to sample such the resultant distribution of values should have specified mean and standard deviation. I can think of rejection sampling to achieve this however that seems to be overkill as rejection sampling tries to match complete distribution while I need only first two moment to be equal to specified value. There have to be easy way to do this and also multiple ways.

The current data have distribution

Min.    1st Qu.     Median     Mean     3rd Qu.    Max.   
0.0020  0.2120      0.3880     0.4298   0.6040     1.0000 

I want it to have mean of 0.2236 and standard deviation of 0.3180463.

Answer 1

Without looking at your data, and without needing to, I can bet 1K points of my reputation that you won't be able to pull a sample that has exactly this mean and exactly this standard deviation. This is just a zero probability event.

Without going into much asking about why you want to do this, here's an approach that will give you an approximate solution. Obtain maximum empirical likelihood weights that take your original data as input, and impose the two constraints for the two moments. Treating these weights as probabilities of selection, you will have a population with the true population mean equal to BLAH, and the true population variance equal to BLAH. You can then either pull the unequal probability samples explicitly, or simply expand your original data by a factor of say $10^4 \times$ {the empirical likelihood weight}, to get a finite population that has approximately the right mean and variance (to about four decimal points).

Answer 2

If you have $n$ sample points $(x_1,\dots,x_n)$ and want a subsample of $m$ points for which the subsample mean and subsample standard deviation approximate some given numbers $\mu^*$ and $\sigma^*$, one possibility is to think about it as an Integer Programming problem. A (possibly nonunique) solution is a vector $b\in\{0,1\}^n$, with $\sum_{i=1}^n b_i=m$, for which each coordinate indicates if the corresponding sample point is or is not included in the subsample, and $b$ minimizes something like $$ \left|\mu^* - \frac{1}{m}\sum_{i=1}^n x_ib_i\right| + \left| \sigma^* - \sqrt{\frac{1}{m-1}\sum_{i=1}^n \left(x_ib_i-\left(\frac{1}{m}\sum_{j=1}^n x_jb_j\right)\right)^2}\right| \, . $$ I don't know if this is feasible for your data set. There are freely available Integer Programming packages.

Answer 3

First of all as many pointed out it impossible to exactly match mean and distribution. Apologies for making not clear in the problem. All I wanted was sampled distribution as close as possible to target mean and SD.

Right now I implemented a rejection sampling which is closely linked to the solution suggested by @Stask. (Let's say) I was targeting a normal distribution $Q(x) = N( \mu, \sigma)$. The empirical distribution of the data $P(x)$ (this can be estimated using kernels or windowing function like Parzen). All need to be done is to draw sample from data with repetition uniform distribution for each sample $y$ if $Q(x)/P(x) > 1$ accept the sample otherwise accept with probability $Q(x)/P(x)$. As more sample will drawn the sample distribution will converge to $Q(x)$.

A layman version of the rejection sampling would be to divide $Q(x)$ in (large number of ) quantiles. For each quantile equal number of samples from the data which are in same range as within the quantile of $Q(x)$ . As number of quantile increases sampled distribution converges to $P(x)$.

Answer 4

A simple approach that may work for others with similar problems is to look at a large number of subsets and save the best matches. You'll need to define what "better" means since you have two objectives (mean & std). In this example, I gave equal weight to variations from either target. You probably would want to play with different subsample sizes to find a large N that still has a close enough mean and std. You'll also want to repeat your analysis with different subsamples to make sure your results are robust.

# Parameters
target_mean = 1
target_std = .5
subsample_size=15
num_subsamples_to_keep = 5
num_subsamples_to_try = 10000

# This scoring metric is kind of arbitrary and could be tuned. 
# For this one, I'm giving equal penalty to deviations from the target mean and the 
# target standard deviation. If fitting the target mean was twice as important as 
# fitting the target std, you could use:
#    return(2*abs(mean-target_mean) + abs(std-target_std))
# Lower scores are better than higher scores. If you redefined the scoring function 
# to make high scores better, then you would need to do a reverse sort in the for loop below. 
calculate_score <- function(mean, target_mean, std, target_std) { 
    return(abs(mean-target_mean) + abs(std-target_std))
}

# Make up some data. This will have a different distribution than what we are looking for.
df = data.frame(x=rnorm(100), y=rnorm(100))

# The matrix saved_subsamples will stores the best subsamples, but it will start off blank.
# A corresponding dataframe below will store the scores (subsample_metrics). 
# We'll keep saved_subsamples and subsample_metrics synchronized so that row 1 in one 
# always corresponds to row 1 of the other. 
saved_subsamples = matrix(nrow=num_subsamples_to_keep, ncol=subsample_size)

# Store the metrics of a subsample's quality in this data frame. 
# We'll keep the rows of this synchronized with saved_subsamples and
# keep them both sorted so the worst saved subsample is always in the last row
# When we find a new subsample that's better than the worst known, we'll save 
# it and resort the matrix and dataframe.
subsample_metrics = data.frame(mean=matrix(NaN, nrow=num_subsamples_to_keep), std=matrix(NaN, nrow=num_subsamples_to_keep), score=matrix(NaN, nrow=num_subsamples_to_keep))

# Try a bunch of different subsamples and keep the best ones
for (attempt in 1:num_subsamples_to_try) {
    # Pick a new subsample and evaluate it. 
    subsample = sample(nrow(df),size=subsample_size)
    subsample_mean = mean(df[subsample,'x'])
    subsample_std = sd(df[subsample,'x'])
    subsample_score = calculate_score(subsample_mean, target_mean, subsample_std, target_std)
    # If we don't have any subsamples yet (worst score is NA), or if this new 
    # subsample is better than our worst saved subsample, then overwrite our 
    # worst subsample from the end and re-order the lists. 
    if( is.na(subsample_metrics[nrow(subsample_metrics),'score']) ||
            subsample_score < subsample_metrics[nrow(subsample_metrics),'score'] ) {
        # Overwrite the worst saved subsample (row 1)
        subsample_metrics[nrow(subsample_metrics),'mean'] = subsample_mean
        subsample_metrics[nrow(subsample_metrics),'std'] = subsample_std
        subsample_metrics[nrow(subsample_metrics),'score'] = subsample_score
        # Transpose the subsample so it fits in one row instead of one column
        saved_subsamples[nrow(saved_subsamples),] = t(subsample)
        # Reorder the saved subsamples (and corresponding metrics) so the worst 
        # scores (the largest) are at the bottom of the list.
        new_order = order(subsample_metrics$score)
        saved_subsamples = saved_subsamples[new_order,]
        subsample_metrics = subsample_metrics[new_order,]
    }
}

# Look at the top subsamples found & their scores
print(subsample_metrics)
print(saved_subsamples)

# Retrieve the top-scoring subsample from the end of the list
top_subsample_indexes = saved_subsamples[1,]
# Use that to subsample the original dataset
top_subsample = df[top_subsample_indexes,]