Compare Raw-BIC and BIC based on RSS.

bic.py

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import norm

def bic_raw(x, k):
    n = len(x)
    return -2 * np.sum(np.log(norm.pdf(x))) + k * np.log(n)

def bic_rss(x, k):
    n = len(x)
    rss = np.sum(x*x)
    return n * np.log(rss/n) + k * np.log(n)
  
n = 100
k = 10

raw_list = []
rss_list = []
for _ in range(1000):
    x = np.random.normal(size=n)
    raw_list.append(bic_raw(x, k))
    rss_list.append(bic_rss(x, k))
raw_mean = np.mean(raw_list)
rss_mean = np.mean(rss_list)
bins = plt.hist(raw_list, bins=50, label="BIC(raw):mean={:.3f}".format(raw_mean))
dbic = n * (1 + np.log(2*np.pi))
plt.hist(rss_list, bins=bins[1]-dbic, label="BIC(rss):mean={:.3f}".format(rss_mean))
plt.legend()
plt.title("Distribution of BIC(raw) & BIC(rss)\nn(1+ln(2π))={:.3f}".format(dbic))
plt.xlabel("BIC")
plt.savefig("dist_bic.png")
plt.close()

d_bic_list = [ra - rs for ra, rs in zip(raw_list, rss_list)]
bins = plt.hist(np.array(d_bic_list), bins=50)
ymax = np.max(bins[0])
dbic_i = 0
for i in range(len(bins[1])-1):
    if bins[1][i] <= dbic and bins[1][i+1] >= dbic:
        dbic_i = i
        break
x = (bins[1][dbic_i] + bins[1][dbic_i + 1]) / 2
width = bins[1][1] - bins[1][0]
plt.bar(x=x, height=bins[0][dbic_i], width=width, color="r", label="Bin of {:.3f}".format(dbic))
plt.legend()
plt.title("Distribution of BIC(raw) - BIC(rss)")
plt.xlabel("BIC(raw) - BIC(rss)")
plt.savefig("dist_dbic.png")