There are lots of different ways of doing it (for a reference, see https://www.mdpi.com/1099-4300/19/11/631 https://www.mdpi.com/1099-4300/19/11/631 ), but this

There are lots of different ways of doing it (for a reference, see https://www.mdpi.com/1099-4300/19/11/631 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 }