Skip to content
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’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Weighted KDE #26

Merged
merged 12 commits into from
Jun 27, 2017
Merged

Weighted KDE #26

Changes from 2 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
ce9ca33
extend tabulate, wrap with weights=ones(data)
axsk Jun 10, 2016
903070f
finish
axsk Jun 10, 2016
3ade534
Merge branch 'master' into HEAD
axsk Oct 13, 2016
0af7754
implement UniformWeights
axsk Oct 13, 2016
fd3a44a
add bivariate weights
axsk Oct 13, 2016
af24d40
change tabulate() signature; use default_weights()
axsk Oct 13, 2016
3ff043e
add length parameter to type
axsk Oct 13, 2016
7ad63b0
add test for weighted univ. kde
axsk May 18, 2017
26529d9
Merge branch 'master' into weights
axsk May 18, 2017
b1491ee
fix test bracket
axsk May 18, 2017
346874f
remove unneeded methods from UniformWeights
axsk May 18, 2017
4a3680a
add bivariate tests
axsk May 18, 2017
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
37 changes: 23 additions & 14 deletions src/univariate.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,29 +67,33 @@ function kde_range(boundary::(@compat Tuple{Real,Real}), npoints::Int)
end

# tabulate data for kde
function tabulate(data::RealVector, midpoints::Range)
ndata = length(data)
function tabulate(data::RealVector, weights::RealVector, midpoints::Range)
npoints = length(midpoints)
s = step(midpoints)

# Set up a grid for discretized data
grid = zeros(Float64, npoints)
ainc = 1.0 / (ndata*s*s)
ainc = 1.0 / (sum(weights)*s*s)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When you switch to WeightsVec, you can use its sum field to avoid recomputing its sum.


# weighted discretization (cf. Jones and Lotwick)
for x in data
for (i,x) in enumerate(data)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Even better: for (x, w) in zip(data, weights).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I totally agree here, but this needs weights to implement the "Iterator protocoll", which is not given for weights::WeightVec or my (WIP) custom UniformWeights type.

k = searchsortedfirst(midpoints,x)
j = k-1
if 1 <= j <= npoints-1
grid[j] += (midpoints[k]-x)*ainc
grid[k] += (x-midpoints[j])*ainc
grid[j] += (midpoints[k]-x)*ainc*weights[i]
grid[k] += (x-midpoints[j])*ainc*weights[i]
end
end

# returns an un-convolved KDE
UnivariateKDE(midpoints, grid)
end

function tabulate(data::RealVector, midpoints::Range)
weights = ones(data)
tabulate(data, weights, midpoints)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could do something like tabulate(data, midpoints; weights=ones(data)), which makes the weighting more apparent (since it's a keyword argument) and removes the need for two separate methods like this.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As you say below, better use the WeightsVec type than a keyword argument for consistency.

end

# convolve raw KDE with kernel
# TODO: use in-place fft
function conv(k::UnivariateKDE, dist::UnivariateDistribution)
Expand Down Expand Up @@ -118,31 +122,36 @@ function conv(k::UnivariateKDE, dist::UnivariateDistribution)
UnivariateKDE(k.x, dens)
end

function uniformweights(data)
n = length(data)
fill(1/n, n)
end

# main kde interface methods
function kde(data::RealVector, midpoints::Range, dist::UnivariateDistribution)
k = tabulate(data, midpoints)
function kde(data::RealVector, weights::RealVector, midpoints::Range, dist::UnivariateDistribution)
k = tabulate(data, weights, midpoints)
conv(k,dist)
end

function kde(data::RealVector, dist::UnivariateDistribution;
boundary::(@compat Tuple{Real,Real})=kde_boundary(data,std(dist)), npoints::Int=2048)
boundary::(@compat Tuple{Real,Real})=kde_boundary(data,std(dist)), npoints::Int=2048, weights=uniformweights(data))

midpoints = kde_range(boundary,npoints)
kde(data,midpoints,dist)
kde(data,weights,midpoints,dist)
end

function kde(data::RealVector, midpoints::Range;
bandwidth=default_bandwidth(data), kernel=Normal)
bandwidth=default_bandwidth(data), kernel=Normal, weights=uniformweights(data))
bandwidth > 0.0 || error("Bandwidth must be positive")
dist = kernel_dist(kernel,bandwidth)
kde(data,midpoints,dist)
kde(data,weights,midpoints,dist)
end

function kde(data::RealVector; bandwidth=default_bandwidth(data), kernel=Normal,
npoints::Int=2048, boundary::(@compat Tuple{Real,Real})=kde_boundary(data,bandwidth))
npoints::Int=2048, boundary::(@compat Tuple{Real,Real})=kde_boundary(data,bandwidth), weights=uniformweights(data))
bandwidth > 0.0 || error("Bandwidth must be positive")
dist = kernel_dist(kernel,bandwidth)
kde(data,dist;boundary=boundary,npoints=npoints)
kde(data,dist;boundary=boundary,npoints=npoints,weights=weights)
end

# Select bandwidth using least-squares cross validation, from:
Expand Down