2
$\begingroup$

I was handed a discrete histogram of a random variable $x$. How do I generate 2,000 continuous samples from the histogram which represent the original random variable $x$?

My first thought is to fit a line to the cdf of the histogram then generate samples under the line using the inverse transform method. Is this a good approach? Do other sensible methods exist?

Clarification

By continuous I mean that the random variable $x$ is continuous, but all I have is a histogram which is a discrete representation of $x$.

Background/Purpose

I have millions of models, each with 40,0000 samples from Bayesian models. I do not have the space to store the 40,0000 samples generated from the these prediction. The goal is to store a function, histogram (the focus of this question), summary statistics, samples, etc. that represent the original distributions. A simple solution would be to store a smaller number of the samples, but I am interested in trying other methods.

$\endgroup$
2
  • $\begingroup$ Could you clarify what you mean by a "continuous sample"? Also, explaining the purpose of this sampling process could help us identify suitable solutions out of the very many possible ones. $\endgroup$
    – whuber
    Commented Dec 22, 2020 at 17:01
  • $\begingroup$ This is getting interesting--but it's turning into a completely different problem. Is your basic concern one of storing a distribution in a compact but reasonably accurate form? If so, how compact and how accurate must it be? What operations must it support efficiently (besides sampling from it)? As an example of one approach that differs from storing histograms or other discrete approximations, see the solutions at stats.stackexchange.com/questions/35220. $\endgroup$
    – whuber
    Commented Dec 22, 2020 at 18:21

2 Answers 2

1
$\begingroup$

Here is an idea that will allow you to capture the 'general shape' of the original large sample using the 512 x-coordinates of the kernel density estimator and its 512 y-coordinates.

Start with $n = 40\,000$ observations from some population, for example a population distributed $\mathsf{Gamma}(\mathrm{shape}=5,\mathrm{rate}=0.01).$

set.seed(2020)
n = 40000
set.seed(2020)
w.orig = rgamma(n, 5, .01)
summary(w.orig);  length(w.orig);  sd(w.orig)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  37.87  338.45  466.78  500.57  626.47 1964.60 
[1] 40000
[1] 222.3676

It's default kernel density estimator (KDE) from R is consists of 512 $(x,y)$ coordinates, which generally capture the shape of the sample more accurately than a histogram. Coordinates summarized below:

density(w.orig)

Call:
        density.default(x = w.orig)

Data: w.orig (40000 obs.);      Bandwidth 'bw' = 23.24

       x                 y            
 Min.   : -31.84   Min.   :4.900e-09  
 1st Qu.: 484.70   1st Qu.:5.862e-06  
 Median :1001.23   Median :9.144e-05  
 Mean   :1001.23   Mean   :4.835e-04  
 3rd Qu.:1517.77   3rd Qu.:8.463e-04  
 Max.   :2034.31   Max.   :1.986e-03  

Then we can take a sample of the x-coordinates in proportion to the y-coordinates to 'restore' the sample. The left-hand panel of the figure below shows a histogram of the original sample along with its KDE. The panel at right shows a histogram of the 'restored' sample.

set.seed(1222)
w.rest = sample(density(w.orig)$x, n, rep=T, p=density(w.orig)$y)
summary(w.rest); length(w.rest);  sd(w.rest)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -23.75  336.10  465.49  498.61  627.22 1969.61 
[1] 40000
[1] 221.6098

par(mfrow=c(1,2))
 hist(w.orig, prob=T, col="skyblue2", main="Original")
  lines(density(w.orig), col="orange")
 hist(w.rest, prob=T, col="cyan", main="Restored")
par(mfrow=c(1,1))

enter image description here

Notes: (1) With its default 'bandwidth', this KDE shows a very few negative plotting points (not noticeable on the plot, perhaps unimportant). If these pose a technical difficulty, you can likely delete either the few $(x,y)$ points with negative $x$'s, or the few negative points in w.rest, with minimal damage to the 'restored' sample.

(2) If storing $512 * 2 = 1024$ coordinate values takes too much space, you could probably get good results by saving every 2nd or 3rd pair of coordinates.

(3) Of course, if you knew that your data are normal, exponential, gamma or from some other family of distributions, then you could estimate the distribution parameters, and use the few sufficient statistics to reconstruct the sample.

(4) Even by saving a subsample of size $1024,$ I don't think you will do nearly as well, trying to use a 'bootstrap' re-sample from that subsample to imitate the sample. (Ordinarily, bootstrapping would use the original full sample.) And if you were to 'save' fewer than 1024 observations. then your imitation would usually be even worse.

set.seed(1234)
w.sub = sample(w.orig, 1024)
w.boot = sample(w.sub, n, rep=T)
summary(w.boot)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   82.5   346.5   482.7   510.1   636.6  1638.9 

par(mfrow=c(1,2))
 hist(w.orig, prob=T, col="skyblue2", main="Original")
 hist(w.boot, prob=T, col="wheat", main="'Bootstrap' Restoration")
par(mfrow=c(1,1))

enter image description here

(5) Finally, your idea to use the histogram might work OK. If you think the histogram of the original data has enough bins so that it represents the original data very well, then here is how to capture histogram information to imitate the original sample. Making a histogram with parameter plot=F gives information about the histogram; the important information is the midpoints of the histogram bars and the bar frequencies or densities.

[You might plot a KDE through the histogram to judge how good the histogram is. Because the histogram above seems to fit well, an imitation based on its information should be useful.]

set.seed(1235)
hist.stats = hist(w.orig, plot=T)
w.hist = sample(hist.stats$mids, n, rep=T, hist.stats$dens)
summary(w.hist)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   50.0   350.0   450.0   502.5   650.0  1950.0 

par(mfrow=c(1,2))
 hist(w.orig, prob=T, col="skyblue2", main="Original")
 hist(w.hist, prob=T, col="wheat", main="Histgraom Restoration")
par(mfrow=c(1,1))

enter image description here

$\endgroup$
0
$\begingroup$

Your approach makes sense. You can fit an empirical distribution using the discrete points you have and then take a random draw from this distribution. This is how bootstrapping works.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.