Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
mutual information pure
function naivehist(x::AbstractVector, ε::Real)
mi, ma = extrema(x)
hs = zeros(Int((ma-mi)÷ε) + 1)
is = zeros(Int, length(x))
for (j, e) in enumerate(x)
xi = Int((e-mi)÷ε + 1)
hs[xi] += 1
is[j] = xi
end
return hs ./ length(x), is
end
function naivehist(x::AbstractVector, y::AbstractVector, ε::Real)
@assert length(x) == length(y)
x₋, x₊ = extrema(x)
y₋, y₊ = extrema(y)
pxy = zeros(Int((x₊ - x₋)÷ε) + 1, Int((y₊ - y₋)÷ε) + 1)
return naivehist!(pxy, x, y, ε)
end
function naivehist!(pxy, x, y, ε)
pxy .= 0.0
x₋ = minimum(x)
y₋ = minimum(y)
for j in 1:length(x)
xv, yv = x[j], y[j]
xi = Int((xv - x₋)÷ε + 1)
yi = Int((yv - y₋)÷ε + 1)
pxy[xi, yi] += 1
end
return pxy ./= length(x)
end
function mutualinfo(x, y, ε; base = MathConstants.e)
px, ix = naivehist(x, ε)
py, iy = naivehist(y, ε)
pxy = naivehist(x, y, ε)
mutualinfo(x, y, px, ix, py, iy, pxy)
end
function mutualinfo(x, y, px, ix, py, iy, pxy, base)
mi = 0.0
x₋, x₊ = extrema(x)
y₋, y₊ = extrema(y)
for (m, xv) in enumerate(x)
xi = ix[m]
for (n, yv) in enumerate(y)
yi = iy[n]
ζ = pxy[xi, yi]
ζ == 0 && continue
mi += ζ * log(base, ζ/(px[xi]*py[yi]))
end
end
return mi
end
using Random: shuffle!
function mutualinfoshuffle(x, y, ε; M = 10_000, base = MathConstants.e)
px, ix = naivehist(x, ε)
py, iy = naivehist(y, ε)
pxy = naivehist(x, y, ε)
return mutualinfoshuffle!(copy(x), copy(y), px, ix, py, iy, pxy, M, base)
end
function mutualinfoshuffle!(x, y, px, ix, py, iy, pxy, M, base)
mi = mutualinfo(x, y, px, ix, py, iy, pxy)
mis = zeros(M)
for k in 1:M
shuffle!(x); shuffle!(y)
naivehist!(pxy, x, y, ε)
mis[k] = mutualinfo(x, y, px, ix, py, iy, pxy, base)
end
return mi, mis
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment