I'm following the third recipe of this answer to estimate the Shannon entropy of my samples using histograms. My expectation was, increasing the sample size should lead to a better estimation of the true entropy. To test that, I sampled from a Gaussian distribution $N(0, \sigma^2)$ with known entropy $H(X) = 0.5[1 + \log(2\pi \sigma^2)]$:

import numpy as np
from scipy.stats import entropy # by default in natural log
import matplotlib.pyplot as plt

sigma = 2
H0 = 0.5*(1 + np.log(2*np.pi*sigma**2)) # theoretical value in nats

Hs = []
ns = np.logspace(2,6,5, dtype=int)
for n in ns:
    X = np.random.normal(0, sigma, size=n)
    nbins = int(n/20)  # <--- This is causes divergence, read UPDATE 2
    hist, bin_edges = np.histogram(X, nbins, density=True)
    Hs.append(entropy(hist)) # will be automatically normalized to sum to 1

plt.plot(ns, Hs, '-x', label='estimation')
plt.hlines(H0, ns[0], ns[-1], linestyle='--', color='k', label='exact H') 
plt.legend()
plt.yticks(np.arange(11))
plt.xscale('log')
plt.xlabel('num samples')
plt.ylabel('sample entropy')

To my surprise, they diverge:

I bet I'm missing something very simple but unfortunately, cannot spot it. Any help is appreciated.

UPDATE

I found a relevant question that is concerned about the case where the probability distribution is either "peaked" or "flat". The upshot of the most voted answer is that one has to normalize the entropy by $\log n$ where $n$ is the number of samples. Although I find this normalization sensible, I still observe an offset between the theoretical value and the estimation:

plt.plot(ns, np.array(Hs)/np.log(ns), '-x', label='normalized estimation') # normalized

What causes this discrepancy?

UPDATE 2

As @sextus-empiricus noted below, scaling the number bins with the sample size leads to a coarse density estimation. In fact, fixing the number of bins, resolves the divergence. However, no matter what number of bins one chooses, there's still some offset between the estimated (normalized or non-normalized) entropy and the theoretical value, as you can see below:

I appreciate it if someone can explain the reason behind this gap as well. The curve displacement by nbins clearly suggest that there's a relationship between the two, but I found the reason elusive.

I'm following the third recipe of this answer to estimate the Shannon entropy of my samples using histograms. My expectation was, increasing the sample size should lead to a better estimation of the true entropy. To test that, I sampled from a Gaussian distribution $N(0, \sigma^2)$ with known entropy $H(X) = 0.5[1 + \log(2\pi \sigma^2)]$:

import numpy as np
from scipy.stats import entropy # by default in natural log
import matplotlib.pyplot as plt

sigma = 2
H0 = 0.5*(1 + np.log(2*np.pi*sigma**2)) # theoretical value in nats

Hs = []
ns = np.logspace(2,6,5, dtype=int)
for n in ns:
    X = np.random.normal(0, sigma, size=n)
    nbins = int(n/20)  # <--- This is causes divergence, read UPDATE 2
    hist, bin_edges = np.histogram(X, nbins, density=True)
    Hs.append(entropy(hist)) # will be automatically normalized to sum to 1

plt.plot(ns, Hs, '-x', label='estimation')
plt.hlines(H0, ns[0], ns[-1], linestyle='--', color='k', label='exact H')
plt.legend()
plt.yticks(np.arange(11))
plt.xscale('log')
plt.xlabel('num samples')
plt.ylabel('sample entropy')

To my surprise, they diverge:

Am I missing something very simple?

I found a relevant question that is concerned about the case where the probability distribution is either "peaked" or "flat". The upshot of the most voted answer is that one has to normalize the entropy by $\log n$ where $n$ is the number of samples. Although I find this normalization sensible, I still observe an offset between the theoretical value and the estimation:

plt.plot(ns, np.array(Hs)/np.log(ns), '-x', label='normalized estimation') # normalized

What causes this discrepancy?

As @sextus-empiricus noted below, scaling the number bins with the sample size leads to a coarse density estimation. In fact, fixing the number of bins, resolves the divergence. However, no matter what number of bins one chooses, there's still some offset between the estimated (normalized or non-normalized) entropy and the theoretical value, as you can see below:

What is the reason behind this gap? The curve displacement by nbins clearly suggest that there's a relationship between the two, but I found the reason elusive.

I'm following the third recipe of this answer to estimate the Shannon entropy of my samples using histograms. My expectation was, increasing the sample size should lead to a better estimation of the true entropy. To test that, I sampled from a Gaussian distribution $N(0, \sigma^2)$ with known entropy $H(X) = 0.5[1 + \log(2\pi \sigma^2)]$:

import numpy as np
from scipy.stats import entropy # by default in natural log
import matplotlib.pyplot as plt

sigma = 2
H0 = 0.5*(1 + np.log(2*np.pi*sigma**2)) # theoretical value in nats

Hs = []
ns = np.logspace(2,6,5, dtype=int)
for n in ns:
    X = np.random.normal(0, sigma, size=n)
    nbins = int(n/20)  # <--- This is causes divergence, read UPDATE 2
    hist, bin_edges = np.histogram(X, nbins, density=True)
    Hs.append(entropy(hist)) # will be automatically normalized to sum to 1

plt.plot(ns, Hs, '-x', label='estimation')
plt.hlines(H0, ns[0], ns[-1], linestyle='--', color='k', label='exact H') 
plt.legend()
plt.yticks(np.arange(11))
plt.xscale('log')
plt.xlabel('num samples')
plt.ylabel('sample entropy')

To my surprise, they diverge:

I bet I'm missing something very simple but unfortunately, cannot spot it. Any help is appreciated.

UPDATE

I found a relevant question that is concerned about the case where the probability distribution is either "peaked" or "flat". The upshot of the most voted answer is that one has to normalize the entropy by $\log n$ where $n$ is the number of samples. Although I find this normalization sensible, I still observe an offset between the theoretical value and the estimation:

plt.plot(ns, np.array(Hs)/np.log(ns), '-x', label='normalized estimation') # normalized

What causes this discrepancy?

UPDATE 2

As @sextus-empiricus noted below, scaling the number bins with the sample size leads to a coarse density estimation. In fact, fixing the number of bins, resolves the divergence. However, no matter what number of bins one chooses, there's still some offset between the estimated (normalized or non-normalized) entropy and the theoretical value, as you can see below:

I appreciate it if someone can explain the reason behind this gap as well.

I'm following the third recipe of this answer to estimate the Shannon entropy of my samples using histograms. My expectation was, increasing the sample size should lead to a better estimation of the true entropy. To test that, I sampled from a Gaussian distribution $N(0, \sigma^2)$ with known entropy $H(X) = 0.5[1 + \log(2\pi \sigma^2)]$:

import numpy as np
from scipy.stats import entropy # by default in natural log
import matplotlib.pyplot as plt

sigma = 2
H0 = 0.5*(1 + np.log(2*np.pi*sigma**2)) # theoretical value in nats

Hs = []
ns = np.logspace(2,6,5, dtype=int)
for n in ns:
    X = np.random.normal(0, sigma, size=n)
    nbins = int(n/20)  # <--- This is causes divergence, read UPDATE 2
    hist, bin_edges = np.histogram(X, nbins, density=True)
    Hs.append(entropy(hist)) # will be automatically normalized to sum to 1

plt.plot(ns, Hs, '-x', label='estimation')
plt.hlines(H0, ns[0], ns[-1], linestyle='--', color='k', label='exact H') 
plt.legend()
plt.yticks(np.arange(11))
plt.xscale('log')
plt.xlabel('num samples')
plt.ylabel('sample entropy')

To my surprise, they diverge:

I bet I'm missing something very simple but unfortunately, cannot spot it. Any help is appreciated.

UPDATE

I found a relevant question that is concerned about the case where the probability distribution is either "peaked" or "flat". The upshot of the most voted answer is that one has to normalize the entropy by $\log n$ where $n$ is the number of samples. Although I find this normalization sensible, I still observe an offset between the theoretical value and the estimation:

plt.plot(ns, np.array(Hs)/np.log(ns), '-x', label='normalized estimation') # normalized

What causes this discrepancy?

UPDATE 2

As @sextus-empiricus noted below, scaling the number bins with the sample size leads to a coarse density estimation. In fact, fixing the number of bins, resolves the divergence. However, no matter what number of bins one chooses, there's still some offset between the estimated (normalized or non-normalized) entropy and the theoretical value, as you can see below:

I appreciate it if someone can explain the reason behind this gap as well. The curve displacement by nbins clearly suggest that there's a relationship between the two, but I found the reason elusive.