# Any good working method for sampling from estimated distribution in R

Suppose I have a number of points $\{x_i\}_{i=1}^N$ and I want to estimate the smooth distribution and sample from it.

I have tried at the moment a lot of things, and all of them show poor performance. I just want to get the function for one and two-dimensional data, for which the histogram looks like the estimated distribution and doesn't give weird results.

Please, can anyone help me to construct a working method, that allows me to create continuous distribution in $\mathbb{R}$ from arbitrary number of points $\{x_i\}_{i=1}^N$ in 1-D and 2-D and then sample from it?

What I did:

I've got good results (long story) with ks package + Metropolis algorithm from mcmc, but unfortunately it crashes my R.

Using density from stats with approxfun to approximate the pdf:

library(MASS)

density_generator<-function(s)
{
# this function returns kernel density estimate built on sample 's'
hpi1<-bw.nrd(s) # calculating h parameter for kernel estimation
# density gives the table of points x and y
# approxfun does linear approximation
fhat.pi1 <- density(x=s,bw=hpi1, n=4096)
fhat.pi1$y[1]=0 # cutting tails fhat.pi1$y[4096]=0
approxfun(fhat.pi1,rule=2) #returning density function
}

## testing the density generator
conditional_density_y<-density_generator(c(-4,1,2,3))

integrate(conditional_density_y, lower = -Inf, upper= Inf) # checking that we're close to 1
# 0.9994542 with absolute error < 8.3e-05

y<-seq(from=-7,to=11, by=0.01)
plot(y,sapply(y,conditional_density_y),pch=".")


Then using Metropolis-Hastings algorithm from mcmc:

library(mcmc)

sample_element<-function(densiy_func,nbatches=1000)
{
# this function samples one value x from the given density function
# using Metropolis - Hastings algorithm.
# Starting point is always 0.5, we're doing 1000 batches,
# after that the histogram should start to look like
# the actual distribution.
metrop(densiy_func,0.5,nbatch=nbatches)$batch[nbatches-1] } out<-replicate(1000,sample_element(conditional_density_y)) hist(out, breaks=100)  The result is poor: Use of brute force with inverse CDF works awful, integration often gives weird results #### Brute force sampling using inverse cdf cdf<-function(x) integrate(conditional_density_y,lower=-Inf,upper=x)$value ## gives 0 on big values
library(GoFKernel)
idf<-inverse(cdf)

out<-replicate(1000,idf(runif(1)))
hist(out, breaks=100)


The better way is to use kernel smooth cdf from another package, it sometimes works, and sometimes gives weird results:

s<-c(1,2,3,4,5)
library(DiagTest3Grp)
library(GoFKernel)
bw<-BW.ref(s, method = "KS-SJ")
cdf2<-function(z){KernelSmoothing.cdf(xx=s, c0=z, bw=bw)} # kernel smoothed cdf
idf2<-inverse(cdf2)
out<-replicate(1000,idf2(runif(1)))
par(mar = rep(2, 4))
hist(out, breaks=100)


If all you want to do is estimate the distrubition of a random sample, $\bf{x}$ = $(x_1, x_2, ..., x_N)$, then just bootstrap (aka nonparametric bootstrapping).
• I need to implement specific Hidden Markov chain algorithm, the chain has continuous state and observation spaces. This algorithm generates samples from initial distribution, gives them weights according to observations and then generates new distribution. The output after $N$ steps is estimated transition and emission distributions. May 21, 2016 at 20:41