# R: estimate joint density and sample from conditional densities

I'm trying to implement an algorithm that estimates time-homogeneous Markov chain with continuous state space.

The probability of transition from state $i$ to state $i+1$ has continuous conditional PDF $\mu(x_{i+1}|x_i)$.

Suppose I have $N$ observations of transitions between states, i.e. $(x_i,x_{i+1}), i=1..N$. I need to obtain continuous PDFs $\mu(x_{i+1}|x_i)$(forward) and $\mu(x_i|x_{i+1})$(backward) to sample from.

What I did: To obtain forward density I have generated joint density $\mu(x_i,x_{i+1})$ from these $N$ observations, and expect to generate marginal densities $f(x_i)$, $g(x_{i+1})$ using first and second values from these pairs. Then my conditional densities should be

$$\mu(x_{i+1}|x_i) = \frac{\mu(x_i,x_{i+1})}{f(x_i)}, \ \ \mu(x_i|x_{i+1}) = \frac{\mu(x_i,x_{i+1})}{g(x_{i+1})}.$$

So the first question is whether it is a correct approach or not?

I've made a test implementation on R using library(ks) that Google gives as the most relevant for joint densities. It works strange, R session crashes if I try to set interval bigger than $(-6,6)$ for $x$.

Thus in my second question I'm asking for some efficient way to do this in $R$, because I suppose to re-evaluate the densities many times during estimation and then somehow sample from them.

# Generating 2-D density sample (from ks tutorial)
library(ks)
set.seed(8192)
samp <- 200
mus <- rbind(c(-2,2), c(0,0), c(2,-2))
Sigmas <- rbind(diag(2), matrix(c(0.8, -0.72, -0.72, 0.8), nrow=2), diag(2))
cwt <- 3/11
props <- c((1-cwt)/2, cwt, (1-cwt)/2)
x <- rmvnorm.mixt(n=samp, mus=mus, Sigmas=Sigmas, props=props)

plot(x)

# Gives unconstrained 2x2 H-matrix for Kernel parametrization
Hpi1 <- Hpi(x=x)

# Kernel density estimate
fhat.pi1 <- kde(x=x, H=Hpi1)

# plot
plot(fhat.pi1)
#-------------------------------------

# Conditional density equals joint density / marginal :
# f(y|x) = f(x,y) / f(x)

h<-hpi(x=x[,1]) # 1-D parameter
f_x<-kde(x=x[,1],h) # 1-D marginal density

# checking that it works
predict(f_x,x=1)
predict(fhat.pi1, x=c(0.5,0.5))

# creating function for conditional density with one fixed coordinate x = 0.5
conditional_density_y <-
function(y){predict(fhat.pi1, x=c(y,0.5))/predict(f_x,x=y)}

# generating range of x
y<-seq(from=-6,to=6, by=0.01) # R crashes if bigger than (-6;6)

# plotting
plot(y,sapply(y,conditional_density_y),pch=".")


UPD: So i tried two ways to sample from it: Metropolis - Hastings algorithm and Samplepdf algorithm

The results are the following:

Here is the actual conditional density (no crash for -7..5)

# creating function for conditional density with one fixed coordinate
conditional_density_y <-
function(y){predict(fhat.pi1, x=c(y,0.5))/predict(f_x,x=y)}

# generating range of x
y<-seq(from=-7,to=5, by=0.01)

# plotting
plot(y,sapply(y,conditional_density_y),pch=".")


# Trying to sample using Metropolis - Hastings algorithm

library(mcmc)

myL <- function(y,x){log(predict(fhat.pi1, x=c(y,x))/predict(f_x,x=y))}

out<-replicate(1000,metrop(myL,0.5,nbatch=100,x=0.5)$batch) hist(out, breaks=100)  The right histogram indicates the interesting bug I observed presence of the initial value 0.5 in the sample, that should be filtered (see below). Also I found that number of transitions of the chain nbatch should not me less than 100. The best thing is to pick 100th observation on each try: # Taking another sample --------------------------------------------------------- # creating function for conditional density with one fixed coordinate conditional_density_y <- function(y){predict(fhat.pi1, x=c(y,1))/predict(f_x,x=y)} # generating range of x y<-seq(from=-7,to=5, by=0.01) # plotting plot(y,sapply(y,conditional_density_y),pch=".") myL <- function(y,x){log(predict(fhat.pi1, x=c(y,x))/predict(f_x,x=y))} # taking only the last value out<-replicate(1000,metrop(myL,1,nbatch=100,x=1)$batch[100])

out2<-out[! out %in% 1] # removing the initial value
hist(out2, breaks=50)


# Trying to sample using Samplepdf algorithm

Works much longer, only possible for sample size = 10, because it's stopping from roundoff errors. Setting subdivisions and accuracy doesn't help.

The question about the initial correctness and optimality of the approach is still actual.