KDE Simulation Exercise

Let’s investigate the properties of histograms and kernel density estimates with simulated data. When simulating data, we know the true distribution that generated the data and can compare the true distribution with our histogram and KDE estimates.

We’ll begin by selecting an interesting distribution to generate our data. To highlight the properties of our density estimates, we’d like something that is multimodal. To acheive this, we will simualte from a mixture of Normal densities. Specifically, suppose our random variable of interest has a probability of \(p = 0.7\) being drawn from a \(N(5,7)\) distribtion (where 17 is the standard deviation, not the variance) and probability \(1-p = 0.3\) of being drawn from a \(N(32,5)\) distribution. Mathematically, we can express this as \[\begin{align} Y & = \mathbb{1}_{\{U < 0.7\}} \, X_1 + (1-\mathbb{1}_{\{U < 0.7\}}) \, X_2 \\ X_1 & \sim N(5,7) \\ X_2 & \sim N(32,5) \\ U & \sim Uniform(0,1). \end{align}\]

We can plot the true distribution of \(Y\) in the following manner.

p = 0.7
mu1 = 5
sig1 = 7
mu2 = 32
sig2 = 5
xLower = min(mu1-3*sig1,mu2-3*sig2)
xUpper = max(mu1+3*sig1,mu2+3*sig2)
xGrid = seq(xLower,xUpper,length=10000)
yDens = p*dnorm(xGrid,mu1,sig1)+(1-p)*dnorm(xGrid,mu2,sig2)
plot(xGrid,yDens,type='l')

We can simulate \(N = \{100,1000,10000,100000\}\) observations from this distribution in the following manner.

NVec = c(100,1000,10000,100000)
YSamp = list()
for(i in 1:length(NVec)){
  ySamp = rep(0,NVec[i])
  x1Samp = rnorm(NVec[i],mu1,sig1)
  x2Samp = rnorm(NVec[i],mu2,sig2)
  uSamp = runif(N)
  ySamp[uSamp < p] = x1Samp[uSamp < p]
  ySamp[uSamp >= p] = x2Samp[uSamp >= p]
  YSamp[[i]] = ySamp
}

Notice that I have written a loop for each value of \(N\), but I avoided writing an inner loop for the simulation of each observation. The vectorized simulation code inside the loop is not necessarily the most efficient method of sampling (think about how you could improve it), but it is far better than writing an inner loop. Note also that YSamp is a list object that stores several vectors (one for each value of \(N\)) of simulations from the distribution above. Since each of the vectors has a different length (length \(N\)), we could not store them as columns in a matrix.

With a sample of \(N\) observations from the distribution in hand, we can compute histograms and kernel density estimates.

nBreaks = 50
par(mfrow=c(2,2))
hist(YSamp[[1]],breaks=nBreaks)
hist(YSamp[[2]],breaks=nBreaks)

hist(YSamp[[3]],breaks=nBreaks)
hist(YSamp[[4]],breaks=nBreaks)

adj = 1
par(mfrow=c(2,2))
plot(density(YSamp[[1]],adjust=adj))
plot(density(YSamp[[2]],adjust=adj))

plot(density(YSamp[[3]],adjust=adj))
plot(density(YSamp[[4]],adjust=adj))

The histograms and KDEs demonstrate that the true density is better estimated when you have more data (higher \(N\)). This is almost always true in estimation: more data is better. If you tinker with the breaks and adjust arguments to the hist and density functions, you can further see how your choice of smoothing can effect the estimates. Evidence of oversmoothing and undersmoothing in both the histograms and KDEs will be especially apparent for the simulations with large \(N\); it will be less apparent for small \(N\) because there is less data to smooth.