import tensorflow as tf

    D_int = tf.cast(D, tf.int32)
    N_int = tf.cast(N, tf.int32)
    if y is None:
        y = silverman_rule_of_thumb(N)

    YDistr = tf.contrib.distributions.MultivariateNormalDiag(loc=tf.zeros(D_int, tf.float32), 
                                                             scale_diag=tf.ones(D_int, tf.float32))
    Y = YDistr.sample(N_int)
    T = 1.0/(2.0*N*tf.sqrt(m.pi*y))

    A0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(X, 0), tf.expand_dims(X, 1)), axis=2)
    A = tf.reduce_sum(phi_sampling(A0/(4*y), D))

    B0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(Y, 0), tf.expand_dims(Y, 1)), axis=2)
    B = tf.reduce_sum(phi_sampling(B0/(4*y), D))

    C0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(X, 0), tf.expand_dims(Y, 1)), axis=2)
    C = tf.reduce_sum(phi_sampling(C0/(4*y), D))

    return T*(A + B - 2*C)
