Chapter 4, Gaussian Entropy

# generate 10000 data points from a gaussian distribution
using Plots, Random, Distributions, StatsBase
using LinearAlgebra, SpecialFunctions
σ = 0.25
d = Normal(0,σ)
x = rand(d, 10000)
histogram(x, fmt = :png, legend=:topright, label="data")

w = fit(Histogram, x).weights

# compute entropy naively 
# can because all weights are non-zero
- sum( (w ./ sum(w)) .* (log.(w ./ sum(w))) )

# closed form expression for entropy in 1D
H(σ) = log( 2 * pi * exp(1) * σ^2) / 2.0

H(σ)

# Ahmed and Gokhale estimate of entropy
# see: http://www.nowozin.net/sebastian/blog/
#       the-entropy-of-a-normal-distribution.html
function entropy_ag(X)
    # X is a (k,n) matrix, samples in columns
    k = size(X,1)
    n = size(X,2)
    C = zeros(k,k)
    for i=1:n
        C += X[:,i]*X[:,i]'
    end
    H = 0.5*k*(1.0 + log(pi)) + 0.5*logdet(C)
    for i=1:k
        H -= 0.5*digamma(0.5*(n+1-i))
    end
    H
end

entropy_ag(x')

# to see how entropy changes with variance
plot(H,0,2, fmt = :png, legend=:bottomright, 
    label="entropy", xlabel="variance")

# generate a positive definite hermitian matrix
n = 2
Σ = rand(n,n)
Σ = Σ + Σ'
Σ = Σ + n * I

d = MvNormal(Σ)

length(d)

mean(d)

cov(d)

# MvNormal's built in estimate of entropy
entropy(d)

# test the Ahmed and Gokhale estimate
x = rand(d, 10000)
entropy_ag(x)

# closed form expression
log( (2 * pi * exp(1))^length(d) * det(cov(d)) ) / 2.0

# A naive estimate of entropy 
# in more than one dimension
function entropy_naive(X,bins=10)
    eps = 0.0000001
    (d,n) = size(X)
    (dim= d, length= n)
    counts = reshape(zeros(Int16, bins^d),
        Tuple([bins for i in 1:d]))
    min = minimum(x)
    max = maximum(x) + eps
    step = (max-min)/(bins)
    for i=1:n
        ind = [Int16(floor((x[j,i]-min)/step))+1 
            for j in 1:d]
        counts[Tuple(ind)...]+=1
    end
    c = counts[counts .> 0]
    ent = -sum( (c ./ sum(c)) .* (log.(c ./ sum(c))) )
    return (entropy = ent, histogram = counts)
end

# naive estimate with carefully selected bin count
entropy_naive(x,14).entropy

heatmap(entropy_naive(x,50).histogram, fmt = :png)

# to see that the naive method is naive
# we try it out with different bin counts
trials = [ entropy_naive(x,bins).entropy for bins in 3:100 ];

plot(trials, fmt = :png, legend=:bottomright, 
    label="entropy", xlabel="bins")