How to simulate random variables according to the law of a pregiven data sample

Question

Say I have the following sample:

-0.38848247;0.21655804;-1.08211969;0.00369104;0.0993393;0.25531731;1.40574444;-1.80271115;-0.58780605;-0.35026458;0.67197532;-1.29654652;-0.58874467;-0.4004176;1.46242829;0.83946654;-0.24861179;-1.604154;1.53555232;1.45589014;-0.52092604;1.65504054;0.12854317;1.66771236;1.47964458;-1.55677722;-0.40070632;-0.41957953;0.37803884;0.45349904;1.18265468;1.43979945;-1.55395031;0.01054229;0.48887945;0.12949194;-0.40283111;1.2998402;-1.02205575;-0.34370088;0.22117962;0.52840463;-0.9584811;1.22482249;-1.51189671;0.372597;1.68446854;-0.74440632;1.33920212;0.18072373;-0.35813474;0.39400846;1.18971633;1.02192759;-1.90422461;0.18459334;0.18096905;-1.73870267;1.57349896;-1.05632536;0.1864611;-1.49696658;0.51070568;-2.25007651;-0.64768552;1.80404194;2.25164576;-0.07925576;-0.83550267;-1.65513631;0.25913869;0.36030077;-0.09006407;-1.64359237;-0.1312756;-0.13790883;-0.09940115;-0.02089164;-0.60924589;-0.05451811;2.11736111;-1.40329353;-0.71788744;-0.45888623;-0.75608368;0.45762458;-0.24299548;-0.29224218;-1.0488731;-0.62028903;-0.08257067;0.92297771;0.2964071;-0.02598973;-0.20439059;-0.25195469;0.20572878;-0.49343988;0.25886695;1.20595313

How can I simulate random variables that would have the same density function as that of the above sample? (i.e same law)

A quick and dirty solution is to perform a uniform random selection out of this sample, but I am interested in more sophisticated methods. Any ides? literature references?

Answer 1

I like the idea you proposed. But if you do want alternatives what about fitting a function and then using rejection sampling based on it? Not that it'd be way different but still.

I'm no expert, but how about something like this? Assume the array arr contains your dataset and op the simulated samples.

Red line is the simulated data.

result<-density(arr,n=256)
nsamp<-3000

data<-c()

data$x<-runif(nsamp,min(arr),max(arr))
	data$y<-runif(nsamp,0,max(result$y))

op<-c()

for(i in seq(1,nsamp))
{
    x<-data$x[i]
	    y<-result$y[which.min(abs(x-result$x))]
	    if( data$y[i]<y)
        {op<-rbind(op,c(x))}
}

    hist(op,20)

Answer 2

I would discourage you from using density estimation in such a small data set. In non-parametric density estimation, the bias is on the same order as the variance, which generally is OK if and only if you have sufficient data that the variance is low. In this case both bias and variance will be high--very slight perturbations of your sample will yield wildly fluctuating estimates of your density, which on average are sufficiently incorrect as to be unuseful.

Instead, you should use a wild bootstrap--that is, instead of using a multinomial weight (sampling with replacement), use a continuous weight, e.g. exp(1). Kosorok's free textbook contains useful information in the bootstrap chapter. This adjustment means that every data point will show up in each bootstrap resample. Each bootstrap resample will also have the same density as the original sample, as you requested.

Answer 3

A first approach (which can be justified asymptotically) could be fitting a nonparametric density estimator (KDE, for example) to your data and then sampling from it. I think it is impossible to sample exactly from the same density, unless it is known.

#Your data
x=c(-0.38848247,0.21655804,-1.08211969,0.00369104,0.0993393,0.25531731,
1.40574444,-1.80271115,-0.58780605,-0.35026458,0.67197532,-1.29654652,
-0.58874467,-0.4004176,1.46242829,0.83946654,-0.24861179,-1.604154,
1.53555232,1.45589014,-0.52092604,1.65504054,0.12854317,1.66771236,
1.47964458,-1.55677722,-0.40070632,-0.41957953,0.37803884,0.45349904,
1.18265468,1.43979945,-1.55395031,0.01054229,0.48887945,0.12949194,
-0.40283111,1.2998402,-1.02205575,-0.34370088,0.22117962,0.52840463,
-0.9584811,1.22482249,-1.51189671,0.372597,1.68446854,-0.74440632,
1.33920212,0.18072373,-0.35813474,0.39400846,1.18971633,1.02192759,
-1.90422461,0.18459334,0.18096905,-1.73870267,1.57349896,-1.05632536,
0.1864611,-1.49696658,0.51070568,-2.25007651,-0.64768552,1.80404194,
2.25164576,-0.07925576,-0.83550267,-1.65513631,0.25913869,0.36030077,
-0.09006407,-1.64359237,-0.1312756,-0.13790883,-0.09940115,-0.02089164,
-0.60924589,-0.05451811,2.11736111,-1.40329353,-0.71788744,-0.45888623,
-0.75608368,0.45762458,-0.24299548,-0.29224218,-1.0488731,-0.62028903,
-0.08257067,0.92297771,0.2964071,-0.02598973,-0.20439059,-0.25195469,
0.20572878,-0.49343988,0.25886695,1.20595313)

# sample size
length(x)

# KDE
plot(density(x))

# Sampling from a KDE
samplekde = function(n,data){
resp = vector()
samp = sample(1:length(data),n,rep=T)
h = density(data)$bw
for(i in 1:n) resp[i] = rnorm(1,data[i],h)
return(resp)
}

# Example
x1 = samplekde(100,x)

plot(density(x1))
points(density(x),col="red",type="l")

# qqplot to check how similar x and x1 are
qqplot(x,x1)

Page 5 of this file

http://www.stat.cmu.edu/~cshalizi/350/lectures/28/lecture-28.pdf

The bandwidth parameter used in my code is the one used by default in the R command density().