0
$\begingroup$

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:

rplot01


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)

rplot2

$\endgroup$
1
$\begingroup$

What's the purpose for doing this?

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).

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Slowpoke May 21 '16 at 20:41
  • $\begingroup$ It is some sort of particle filter with its sampling-importance-resampling algorithm. At the moment I'm studying the very first problem of sampling and generating new distribution. $\endgroup$ – Slowpoke May 21 '16 at 20:45
  • $\begingroup$ The original algorithm uses tree density estimation that I can't implement, but authors mention that kernel density estimation is also possible. However, I hadn't figured out yet how to add weights to elements during resampling in the case of kernel estimation. $\endgroup$ – Slowpoke May 21 '16 at 20:48
  • $\begingroup$ As far as I understood, in bootstrapping each new sample is a resample of original data. In this algorithm I'm generating some (continuous) density estimate each time and pick new samples from it. Can you please explain, how I should do bootstrap in my case? $\endgroup$ – Slowpoke May 21 '16 at 21:12
  • $\begingroup$ I'm not very familiar with Hidden Markov chains, so I can't really offer much advice. I think your post is missing some context as to what you are trying to accomplish, that's why I suggested bootstrapping if all you wanted to do is estimate a distribution. However, if you're trying to sample new samples, then you'll have to use MC. However, I'm curious as to how your original data was generated. If there is a distribution to that, then you can attach a reference prior and find the posterior sample your way out of the density estimation problem. $\endgroup$ – Jon May 21 '16 at 21:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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