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)
