Multivariate non-parametric density estimation with many missing values Apologies in advance if any of my terminology here is wrong, I'm not an expert in statistics. If I've made any mistakes, let me know and I'll correct them.
The task I'm looking for some advice on approaching is parameterless density estimation (particularly kernel density estimation, but if somebody wants to suggest an entirely different approach I'd be interested in that too) in which there are many (100s or 1000s) of variables, but for each datapoint, I only know a small number of them (perhaps less than 100). 
The known variables are not flatly distributed, and in fact the value of one variable (known or unknown) may have an effect on how likely another variable is to be observed. Whether or not a particular variable is observed may also have an effect on how likely another variable is to be observed.
Assume here that I have an extremely large number of datapoints. It should also be mentioned that the more often a variable, or combination of variables, is observed, the more important it is for me to be able to predict it. So, for example, if there's a certain value of variable X for which I never observe variable Y, I don't care about being able to estimate Y for that value of X.
Bandwidth selection in this situation is also something I could use some guidance on, for the kernel approach. I'd be happy to hear about either theoretical approaches or R packages.
 A: In order to fit a multivariate kernel density estimate (KDE) with missing data (as it sounds like you are) you need to impute the missing data. Alternatively, you could try using Naive Bayes, which would allow you to fit a unique kernel density estimate per feature, per class, which would alleviate your missing data issue and not require any imputation. Note, Naive Bayes rests on the assumption that your features are independent (although often times we may make this assumption even when it isn't so and the method still works well). 
For the first approach, fitting a multivariate KDE with missing data, there are many imputation approaches available, some of which were designed for this purpose. This well-cited 2014 paper compares a few different approaches: Influence of Missing Values Substitutes on Multivariate Analysis of Metabolomics Data, Gromski et al. 
For the Naive Bayes approach, you'll need to make the assumption that your feautures are independent, thus redefining the original Bayes Rule:
\begin{equation}
P_M(t_k|\textbf{x}) = \frac{p(\textbf{x}|t_k) }{p(\textbf{x})} p(t_k)
\label{eq:nbmultimcprob}
\end{equation}
where $\textbf{x}={x_1, x_2,...,x_N}$ is your vector of features, $t_k$ is class $k$ (and implies a positive class instance $t_k=1$), $p(\textbf{x}|t_k)$ is the likelihood, which you can think of as the multivariate KDE across all features, $p(t_k)$ is the class prior (often based on class frequency of the training set), and  $p(\textbf{x})=\sum_{k=1}^K p\left(t_{k}\right) \prod_{i=1}^{N} p(x_i|t_k)$ is the evidence for $K$ transient classes. Feature independence converts this to Naive Bayes which is defined as:
\begin{equation}
    P_M(t_k|x_1, x_2,..., x_N) = \frac{ \prod_{i=1}^{N} p(x_i|t_k) } {p(\textbf{x})} p(t_k)
\end{equation}
where $p(\textbf{x}|t_k)$  has been replace by the product of probability densities from each feature's KDE: $\prod_{i=1}^{N} p(x_i|t_k)$. 
Lastly, you mention bandwidth-setting in your question. As is mentioned in this 2012 paper: Classification Using Kernel Density Estimates, Ghosh et al., bandwidth is often set uniquely for each class (and in the Naive Bayes setting, perhaps for each feature and each class). Bandwidths are often optimized using a grid search over a range of valid values (which is dictated by the values of the dataset, and whether or not you scale the data). 
