There are lots of different ways of doing it (for a reference, see https://www.mdpi.com/1099-4300/19/11/631 ), but this was the basic approach I was trying that I gave up on because of performance concerns:

fn fast_kde(data: &[f64], grid_size: usize, bandwidth_opt: Option<f64>) -> Vec<f64> {
let n = data.len() as f64;
let min = *data.iter().min_by(|x, y| x.partial_cmp(y).unwrap()).unwrap();
let max = *data.iter().max_by(|x, y| x.partial_cmp(y).unwrap()).unwrap();
let std_dev = (data.iter().map(|&x| x.powi(2)).sum::<f64>() / n - (data.iter().sum::<f64>() / n).powi(2)).sqrt();
// Optimal bandwidth using Silverman's rule or use provided bandwidth
let bandwidth = bandwidth_opt.unwrap_or(1.06 * std_dev * n.powf(-0.2));
// Reflecting data at boundaries
let range = max - min;
let mut reflected_data = data.to_vec();
reflected_data.extend(data.iter().map(|&x| 2.0 * min - x));
reflected_data.extend(data.iter().map(|&x| 2.0 * max - x));
let fft = Radix4::new(grid_size, FftDirection::Forward); // Use the enum
// Create a grid and kernel
let mut kernel = vec![Complex::new(0.0, 0.0); grid_size];
for i in 0..grid_size {
let x = (i as f64 / grid_size as f64) * 2.0 * std::f64::consts::PI;
kernel[i] = Complex::new((-0.5 * (x / bandwidth).powi(2)).exp(), 0.0);
}
// Compute FFT of the kernel
fft.process(&mut kernel);
// Compute the FFT of the data (with kernel applied)
let mut data_fft = vec![Complex::new(0.0, 0.0); grid_size];
for value in reflected_data {
let index = ((value - min) / range * grid_size as f64).floor() as usize % grid_size;
data_fft[index].re += 1.0;
}
fft.process(&mut data_fft);
// Convolve in frequency space
let convolved: Vec<Complex<f64>> = kernel.iter().zip(data_fft.iter()).map(|(k, d)| k * d).collect();
// Inverse FFT
let mut inverse_fft = convolved;
let ifft = Radix4::new(grid_size, FftDirection::Inverse); // Use the enum
ifft.process(&mut inverse_fft);
// Normalize
let sum: f64 = inverse_fft.iter().map(|c| c.re).sum();
inverse_fft.iter().map(|c| c.re / sum).collect()
}

fn entropy(data: &Array1<f64>, bandwidth: f64) -> f64 {
assert!(bandwidth > 0.0, "Bandwidth must be greater than 0");
let n = data.len() as f64;
let normalizing_constant = 0.5 * (2.0 * std::f64::consts::PI * std::f64::consts::E * bandwidth * bandwidth).ln();
let entropy_estimate: f64 = data.iter().map(|&x| {
let p_x = kernel_density(data.as_slice().unwrap(), x, bandwidth);
if p_x > 0.0 {
-p_x * (p_x.ln())
} else {
0.0
}
}).sum();
// Divide by the number of data points and add the normalizing constant
(entropy_estimate / n) + normalizing_constant
}

fn mutual_information_kde(x: &Array1<f64>, y: &Array1<f64>, bandwidth: f64) -> f64 {
assert_eq!(x.len(), y.len());
let n = x.len() as f64;
let x_slice = x.as_slice().unwrap();
let y_slice = y.as_slice().unwrap();
let mi: f64 = x_slice.par_iter().enumerate().map(|(i, &xi)| {
let p_x = kernel_density(x_slice, xi, bandwidth);
let p_y = kernel_density(y_slice, y_slice[i], bandwidth);
let p_xy = joint_kernel_density(x_slice, y_slice, xi, y_slice[i], bandwidth);
if p_xy > 0.0 {
(p_xy / (p_x * p_y)).ln() / n
} else {
0.0
}
}).sum();
let h_x = entropy(x, bandwidth);
let h_y = entropy(y, bandwidth);
let nmi = 2.0 * mi / (h_x + h_y);
nmi
}