This is the code I create for an estimation. But I don’t know what is the command to select bandwidth automatically instead of putting value manually.
<code>import numpy as np
from scipy.stats import norm
from scipy.integrate import quad
X1 = np.array([12, 10, 3, 3, 7, 3, 2, 28, 2, 3, 12, 9, 16, 3, 9, 3, 2, 5, 2])
X2 = np.array([16, 15, 16, 9, 10, 15, 26, 30, 17, 6, 15, 17, 19, 6, 11, 15, 15, 4, 8])
h1n =0.8713
t1 =0
t2 =0
def gaussian_kernel(u):
return norm.pdf(u)
integrals = np.array([h1n*quad(gaussian_kernel, (t2-Yi)/h1n, np.inf)[0] for Yi in X2])
def integrand(x1, t2, h1n):
kernel_values = np.array([gaussian_kernel((x1 - X1k) / h1n) for X1k in X1])
prod_kernel=kernel_values*integrals
sum_kernel = np.sum(prod_kernel)
return (sum_kernel)**2
result, _ = quad(lambda x1: integrand(x1, t2, h1n), t1, np.inf)
result = -0.5 * result
print("The result of the integral is:", result)
</code>
<code>import numpy as np
from scipy.stats import norm
from scipy.integrate import quad
X1 = np.array([12, 10, 3, 3, 7, 3, 2, 28, 2, 3, 12, 9, 16, 3, 9, 3, 2, 5, 2])
X2 = np.array([16, 15, 16, 9, 10, 15, 26, 30, 17, 6, 15, 17, 19, 6, 11, 15, 15, 4, 8])
h1n =0.8713
t1 =0
t2 =0
def gaussian_kernel(u):
return norm.pdf(u)
integrals = np.array([h1n*quad(gaussian_kernel, (t2-Yi)/h1n, np.inf)[0] for Yi in X2])
def integrand(x1, t2, h1n):
kernel_values = np.array([gaussian_kernel((x1 - X1k) / h1n) for X1k in X1])
prod_kernel=kernel_values*integrals
sum_kernel = np.sum(prod_kernel)
return (sum_kernel)**2
result, _ = quad(lambda x1: integrand(x1, t2, h1n), t1, np.inf)
result = -0.5 * result
print("The result of the integral is:", result)
</code>
import numpy as np
from scipy.stats import norm
from scipy.integrate import quad
X1 = np.array([12, 10, 3, 3, 7, 3, 2, 28, 2, 3, 12, 9, 16, 3, 9, 3, 2, 5, 2])
X2 = np.array([16, 15, 16, 9, 10, 15, 26, 30, 17, 6, 15, 17, 19, 6, 11, 15, 15, 4, 8])
h1n =0.8713
t1 =0
t2 =0
def gaussian_kernel(u):
return norm.pdf(u)
integrals = np.array([h1n*quad(gaussian_kernel, (t2-Yi)/h1n, np.inf)[0] for Yi in X2])
def integrand(x1, t2, h1n):
kernel_values = np.array([gaussian_kernel((x1 - X1k) / h1n) for X1k in X1])
prod_kernel=kernel_values*integrals
sum_kernel = np.sum(prod_kernel)
return (sum_kernel)**2
result, _ = quad(lambda x1: integrand(x1, t2, h1n), t1, np.inf)
result = -0.5 * result
print("The result of the integral is:", result)
I want to use Sheather Jones bandwidth selection in my proposed estimation.