Skip to content

gh-115532: Minor tweaks to kde() #117897

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We鈥檒l occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rhettinger merged 3 commits into from
Apr 15, 2024
Merged

gh-115532: Minor tweaks to kde() #117897

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8581eb0
Change variable name to match the supporting reference.
rhettinger Apr 13, 2024
04b4f25
Auto detect size changes
rhettinger Apr 6, 2024
3fe2c65
Add kernel invcdf for the parabolic kernel
rhettinger Apr 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion Doc/library/statistics.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1163,7 +1163,7 @@ accurately approximated inverse cumulative distribution function.
.. testcode::

from random import choice, random, seed
from math import sqrt, log, pi, tan, asin
from math import sqrt, log, pi, tan, asin, cos, acos
from statistics import NormalDist

kernel_invcdfs = {
Expand All @@ -1172,6 +1172,7 @@ accurately approximated inverse cumulative distribution function.
'sigmoid': lambda p: log(tan(p * pi/2)),
'rectangular': lambda p: 2*p - 1,
'triangular': lambda p: sqrt(2*p) - 1 if p < 0.5 else 1 - sqrt(2 - 2*p),
'parabolic': lambda p: 2 * cos((acos(2*p-1) + pi) / 3),
'cosine': lambda p: 2*asin(2*p - 1)/pi,
}

Expand Down
34 changes: 23 additions & 11 deletions Lib/statistics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -919,53 +919,53 @@ def kde(data, h, kernel='normal', *, cumulative=False):
sqrt2pi = sqrt(2 * pi)
sqrt2 = sqrt(2)
K = lambda t: exp(-1/2 * t * t) / sqrt2pi
I = lambda t: 1/2 * (1.0 + erf(t / sqrt2))
W = lambda t: 1/2 * (1.0 + erf(t / sqrt2))
support = None

case 'logistic':
# 1.0 / (exp(t) + 2.0 + exp(-t))
K = lambda t: 1/2 / (1.0 + cosh(t))
I = lambda t: 1.0 - 1.0 / (exp(t) + 1.0)
W = lambda t: 1.0 - 1.0 / (exp(t) + 1.0)
support = None

case 'sigmoid':
# (2/pi) / (exp(t) + exp(-t))
c1 = 1 / pi
c2 = 2 / pi
K = lambda t: c1 / cosh(t)
I = lambda t: c2 * atan(exp(t))
W = lambda t: c2 * atan(exp(t))
support = None

case 'rectangular' | 'uniform':
K = lambda t: 1/2
I = lambda t: 1/2 * t + 1/2
W = lambda t: 1/2 * t + 1/2
support = 1.0

case 'triangular':
K = lambda t: 1.0 - abs(t)
I = lambda t: t*t * (1/2 if t < 0.0 else -1/2) + t + 1/2
W = lambda t: t*t * (1/2 if t < 0.0 else -1/2) + t + 1/2
support = 1.0

case 'parabolic' | 'epanechnikov':
K = lambda t: 3/4 * (1.0 - t * t)
I = lambda t: -1/4 * t**3 + 3/4 * t + 1/2
W = lambda t: -1/4 * t**3 + 3/4 * t + 1/2
support = 1.0

case 'quartic' | 'biweight':
K = lambda t: 15/16 * (1.0 - t * t) ** 2
I = lambda t: 3/16 * t**5 - 5/8 * t**3 + 15/16 * t + 1/2
W = lambda t: 3/16 * t**5 - 5/8 * t**3 + 15/16 * t + 1/2
support = 1.0

case 'triweight':
K = lambda t: 35/32 * (1.0 - t * t) ** 3
I = lambda t: 35/32 * (-1/7*t**7 + 3/5*t**5 - t**3 + t) + 1/2
W = lambda t: 35/32 * (-1/7*t**7 + 3/5*t**5 - t**3 + t) + 1/2
support = 1.0

case 'cosine':
c1 = pi / 4
c2 = pi / 2
K = lambda t: c1 * cos(c2 * t)
I = lambda t: 1/2 * sin(c2 * t) + 1/2
W = lambda t: 1/2 * sin(c2 * t) + 1/2
support = 1.0

case _:
Expand All @@ -974,27 +974,39 @@ def kde(data, h, kernel='normal', *, cumulative=False):
if support is None:

def pdf(x):
n = len(data)
return sum(K((x - x_i) / h) for x_i in data) / (n * h)

def cdf(x):
return sum(I((x - x_i) / h) for x_i in data) / n

n = len(data)
return sum(W((x - x_i) / h) for x_i in data) / n


else:

sample = sorted(data)
bandwidth = h * support

def pdf(x):
nonlocal n, sample
if len(data) != n:
sample = sorted(data)
n = len(data)
i = bisect_left(sample, x - bandwidth)
j = bisect_right(sample, x + bandwidth)
supported = sample[i : j]
return sum(K((x - x_i) / h) for x_i in supported) / (n * h)

def cdf(x):
nonlocal n, sample
if len(data) != n:
sample = sorted(data)
n = len(data)
i = bisect_left(sample, x - bandwidth)
j = bisect_right(sample, x + bandwidth)
supported = sample[i : j]
return sum((I((x - x_i) / h) for x_i in supported), i) / n
return sum((W((x - x_i) / h) for x_i in supported), i) / n

if cumulative:
cdf.__doc__ = f'CDF estimate with {h=!r} and {kernel=!r}'
Expand Down