Hello Everyone,

I would like to reproduce an algorithm that using Bayesian methods to select the bandwidth of a local constant estimator under the context of the nonparametric regression. The basic idea is simply using a multivariate kernel which has the following form:

Where x is a p dimensional vector and H is a bandwidth matrix and for simplicity, just assume it is a diagonal matrix with dimensions of p*p.

The leave-one-out NW estimator is simply:

If we assume y_i \sim Normal (\hat {m}_{-i}(x_i), \sigma ^2) then we can actually write down the likelihood, and if we treat the h_{1,...,p} as unknown paramters we can actually get their point estimates by maximizing the likelihood. Now I just want to implement this in a Bayesian framework so I just write the following stan code:

```
data{
int n_sample; // number of patients
int p; // number of dimensions
matrix[n_sample, p] X; // covaraite matrix
real Y[n_sample]; // vector of responses
real Yf[n_sample]; // same as y but for the loop
}
parameters {
vector [p] bandwidth; // bandwidth
real ysd; // the sd for y , assuming the sd is the same bewtween each treatment, i.e. homogeneity of varaince
}
transformed parameters {
matrix[p, p] D;
// latent continuous response
D = diag_matrix(bandwidth);
}
model {
real m_hat[n_sample]; // latent continuous responses
matrix[n_sample,n_sample] density;
matrix[n_sample, n_sample] reg;
real denominator[n_sample];
real numerator[n_sample];
for (j in 1:p) bandwidth[j]~ uniform(0.0001, 10000 ); // prior on effect sizes of A and B
1/ysd ~ inv_gamma(1, 1); // prior for the variance
for (i in 1:n_sample) {
for (j in 1: n_sample) {
density[i,j] = (2 * 3.14)^(-p/2) * sqrt(1/prod(bandwidth))
* exp( (-1/2)
* ( (X[i,]-X[j,]) * inverse(D) * (X[i,]-X[j,])' ) );
reg[i,j] = density[i,j] * Yf[j] ;
}
denominator[i]=sum(density[i,]) -density[i,i];
numerator[i] = sum(reg[i,])-reg[i,i];
m_hat[i]= (numerator[i]) / (denominator[i]);
Y[i] ~ normal (m_hat[i], ysd);
}
}
```

However, the issue is that I keep getting sampling rejections and the chain just get stuck there for hoursâŚ I think I get something wrong in the loop in the âmodelâ part where I need to code the leave one out estimator. Can anyone give me some ideas about where I did wrong on this? Thanks!!