# How to sample from a distribution so that mean of samples equals expected value?

Given a random variable $X$, how do I have obtain $N$ random variates of $X$ so that the mean value of my samples equals the expected value of $X$?

E.g. let $X$ have uniform distribution on the interval $[0, 1]$, and draw $N = 10000$ samples. We can do this in R using:

set.seed(5)
tmp_vec <- runif(10000, 0, 1)


However, mean(tmp_vec) returns 0.5018471. I would like it so that the mean is 0.5, matching the expected value. Should I rescale the sampled values of tmp_vec?

• I have a function $\sum(f(x))$ over a vector $x$. I want to make a statement concerning how this sum changes when $x$ is randomly perturbed, so that the perturbations follow a known distribution. To do this I am going to sample the sum repeatedly over different draws from the distribution, however, I don't want to have to factor in the variance of the draws into the computation of the sum. – Alex Nov 8 '16 at 1:44
• Or rather, I might sample the distribution once, and then permute the pertubations over the vector $x$. – Alex Nov 8 '16 at 1:50

## 3 Answers

On a related question, namely how to simulate an iid random sample $(X_1,\ldots,X_n)\sim f$, under the constraint that its sum $X_1+\cdots+X_n$ is fixed to an arbitrary value $s_0$, I showed that the density of that sample is given by $$(X_1,\ldots,X_n)\sim f(x_1)\cdots f(x_{n-1})f(s_0-x_1-\cdots-x_n)\mathbb{I}_{s_0-x_1-\cdots-x_{n-1}}(x_n)$$

Note: Under the constraint the sample is no longer independent but the observations are identically distributed, despite the apparent asymmetry in the above.

For instance, a Uniform sample with fixed average $0.5$ would have the joint density $$(U_1,\ldots,U_N)\sim \prod_{i=1}^{N-1}\mathbb{I}_{(0,1)}(u_i)\mathbb{I}_{1/2}(u_1+\cdots+u_N)$$ which can simulated by Gibbs sampling or another MCMC algorithm like RWMH. When using Gibbs sampling, the $N-1$ first coordinates of the sample can be simulated one at a time with $$U_i|U_{-i}\sim \mathbb{I}_{(0,1)}(u_i) \mathbb{I}_{(0,1)}(s_0-u_1-\cdots-u_i-\cdots-u_{n-1})=\mathbb{I}_{(\max\{0,s_0-1-\sum_{j\ne i,n}u_j\},\min\{1,s_0-\sum_{j\ne i,n}u_j\})}(u_i)$$ An R implementation looks like this code:

n=3;T=1e4
s0=.5 #fixed average
sampl=matrix(s0,T,n)
for (t in 2:T){
sampl[t,]=sampl[t-1,]
for (i in 1:(n-1)){
sampl[t,i]=runif(1,
min=max(0,n*s0-sum(sampl[t,c(-i,-n)])-1),
max=min(1,n*s0-sum(sampl[t,c(-i,-n)])))
sampl[t,n]=n*s0-sum(sampl[t,-n])}}


with the following marginals on the $U_i$'s:

Note: One can easily modify the above R code to impose an average of $s_0=0.05$ or $s_0=0.975$ on the sample.

• Thanks! Why is this different from just taking the sample, and then rescaling it? – Alex Nov 8 '16 at 9:44
• If you rescale your sample, it no longer has the original distribution. Nor does it enjoy the proper conditional distribution. You thus end up with an usually unspecified distribution. Take e.g. a Normal sample. When rescaled to have fixed average and empirical variance, it is no longer Normal. – Xi'an Nov 8 '16 at 10:05

I am not sure that you can achieve that. runif is a pseudo-random number generator that produces uniformly distributed random variables. It does not produce random variable with an exact mu that is equal to what you specified. 0.0518 is a pretty good mu. The same is true for the family of functions including rnorm, rexp, rpois etc.

but try this, it should work:

runif2<-function(n,min,max){((min+max)/2)+sqrt((max-min)/12)*scale(runif(n))}
r <- runif2(100,0,1)
mean(r) #0.5


Just keep in mind that this sample is now different than runif because you used 'scale'. It does not have the same properties as your original distribution and it may or may not pass tests for randomness.

Only an infinite sized sample would guarantee an exact result and that is not possible. However, you can set a limit on how close to the "real" mean you would like to get or which is the smallest deviation you are ready to accept in 95% of the cases. The key to that ist the standard error of the mean.

$SE = \frac{s}{sqrt(n)}$

with s being the standard deviation (for large $n$ there is not much point in distinguishing between sample of distribution) you can resolve that to

$n=(\frac{s}{SE})^2$

and thus compute the minimum number fo any given Standard Error that you are willing to accept.

The Wikipedia is actually quite good at this entry: https://en.wikipedia.org/wiki/Standard_error