# Normal distribution
Gaussian = lambda a, b, c : lambda x : a * exp(((x - b) / c)**2 / -2);
Normal = lambda mean, stdev : Gaussian(1 / (sqrt(2 * pi) * stdev), mean, stdev);
Normal_Variance = lambda mean, variance : Normal(mean, sqrt(stdev))
Normal_Standard = Normal(0, 1);
|
|
# Kernel density estimation
KernelDensity = lambda kernel, bandwidth, data : lambda x : sum(map(lambda element : kernel((x - element)**2 / bandwidth), data)) / (len(data) * bandwidth)
KernelDensity_StandardNormalKernel = lambda bandwidth, data : KernelDensity(Normal_Standard, bandwidth, data)
|
|
# Example 1 - As found on Wikipedia at http://en.wikipedia.org/wiki/Kernel_density_estimation
data = [-2.1, -1.3, -0.4, 1.9, 5.1, 6.2];
mean = sum(data) / len(data);
stdev = sqrt(sum(map(lambda element : (element - mean)**2, data)) / len(data));
p = point(map(lambda element : [element, 0], data));
kd1 = plot(KernelDensity(Normal_Variance(0, 2.25), 0.6, data), [min(data) - stdev, max(data) + stdev], color='red') + text('h=0.6', [1.1, 0.16], color='red');
kd2 = plot(KernelDensity(Normal_Variance(0, 2.25), 0.9, data), [min(data) - stdev, max(data) + stdev], color='green') + text('h=0.9', [1.1, 0.15], color='green');
kd3 = plot(KernelDensity(Normal_Variance(0, 2.25), 1.2, data), [min(data) - stdev, max(data) + stdev], color='blue') + text('h=1.2', [1.1, 0.14], color='blue');
show(p + kd1 + kd2 + kd3);
|
|
# Example 2
data = [-2.9, -2.5, -1.5, 2.0, 3.0, 3.2, 3.8, 3.9, 4.5];
mean = sum(data) / len(data);
stdev = sqrt(sum(map(lambda element : (element - mean)**2, data)) / len(data));
p = point(map(lambda element : [element, 0], data));
kd1 = plot(KernelDensity_StandardNormalKernel(0.5, data), [min(data) - stdev, max(data) + stdev], color='red') + text('h=0.5', [-2.0, 0.34], color='red');
kd2 = plot(KernelDensity_StandardNormalKernel(1.3, data), [min(data) - stdev, max(data) + stdev], color='green') + text('h=1.3', [-2.0, 0.32], color='green');
kd3 = plot(KernelDensity_StandardNormalKernel(3.0, data), [min(data) - stdev, max(data) + stdev], color='blue') + text('h=3.0', [-2.0, 0.30], color='blue');
show(p + kd1 + kd2 + kd3);
|
|
