Estimation of the density/distribution Let $(x_i,y_i,z_i)_{i=1,\dots,n}$ be an i.i.d. sample of $(X,Y,Z)$. How one can estimate the following object 
$$\int_{-\infty}^xf(\bar x,y|z)\mathrm{d}\bar x$$
where $f(x,y|z)$ is a density of $X,Y$ conditional on $Z$.
We know that $F(x) = \int_{-\infty}^xf(\bar x)\mathrm{d}x$ can be estimated using $\hat F(x) = \frac{1}{n}\sum_{i=1}^n1[x_i\leq x]$. We also know that $\hat f(x,y|z) = \frac{\frac{1}{h_{xn}h_{yn}}\sum_{i=1}^nK\left(\frac{x-x_i}{h_{xn}}\right)K\left(\frac{y-y_i}{h_{yn}}\right)K\left(\frac{z-z_i}{h_{zn}}\right)}{\sum_{i=1}^nK\left(\frac{z-z_i}{h_{zn}}\right)}$
Can I construct the estimator for the above mentioned object like this?
$$\frac{\frac{1}{h_{yn}}\sum_{i=1}^n1[x_i\leq x]K\left(\frac{y-y_i}{h_{yn}}\right)K\left(\frac{z-z_i}{h_{zn}}\right)}{\sum_{i=1}^nK\left(\frac{z-z_i}{h_{zn}}\right)}$$
where $K$ is a kernel function $1$ is an indicator variable, $h$ is a bandwidth and $n$ is a sample size.
 A: I guess you have a sample from $(X,Y,Z)$. By definition, we have that
$$f(x,y\vert z)=\dfrac{f(x,y,z)}{f(z)}.$$
Using the joint sample from $(X,Y,Z)$, construct a nonparametric estimator (e.g. KDE) of their corresponding density $\varphi_{X,Y,Z}$. Using the marginal sample from $Z$ construct a nonparametric (KDE) estimator of the density of $Z$, $\varphi_Z$. Plug these estimators in the in integral 
$$\int_{-\infty}^xf(\bar x,y|z)\mathrm{d}\bar x \approx \int_{-\infty}^x \dfrac{\varphi_{X,Y,Z}(\bar{x},y,z)}{\varphi_Z(z)}\mathrm{d}\bar x,$$
and voilà. The integral can be calculated using quadrature or Monte Carlo methods.
Note that the bandwidth parameter might not be the same for all the entries as you are suggesting. In addition, possible dependencies of $X,Y,$ and $Z$ must be taken into account. 
A: Bayes, 
Actually, you don't need to estimate separate KDEs for f(x,y,z) and f(z) to get f(x,y,z). You can just estimate the joint $f(x,y,z)$ by KDE that uses the Gaussain kernels and then analytically compute $f(x,y|z)$ for any selected value of z. Then, as Pigeon already explained, just integrate over the x. But, if you are using the Gaussian kernels, you can do this part analytically as well. So you are in effect completely avoiding Monte Carlo. You will, however, require samples from the joint pdf f(x,y,z) and I assume you can get those. Then you'll need a good multivariate KDE that uses Gaussian kernels. Perhaps you can try this one: http://www.mathworks.com/matlabcentral/fileexchange/41187-fast-kernel-density-estimator-multivariate
Note that your equation for $\hat f(x,y|z)$ assumes independance of x,y,z, but I don't think that your initial definition makes that assumption. Anyway, if you're able to get the i.i.d. samples from the joint $f(x,y,z)$, then the above KDE will take into acount the possible correlation accross x,y and z .
You can find the neccesary equations on how to compute conditionals and marginals of multivariate Gaussians in textbooks or check the net.
