Utilities.py 文件源码

def bic_gaussian_kernel(chain_fname,data_length):
    """
    Bayesian information criterion
    Only valid if data_length >> n_params

    # Note that bandwidth of kernel is set to 1 universally
    """
    chain = np.load(chain_fname + '.npy')
    n_params = np.shape(chain)[-1]
    samples = chain.reshape((-1,n_params))
    kde = KernelDensity(kernel='gaussian',bandwidth=1).fit(samples)

    # Best fit = medians of the distribution
    medians = np.median(samples,axis=0) # shape = (n_params,)
    medians = medians.reshape(1,-1)     # Reshape to (1,n_params)

    log10_L = float(kde.score_samples(medians)) 
    lnL = log10_L /np.log10(2.7182818)

    return -2*lnL + n_params*np.log(data_length)


###############################################################################
# Process chain
###############################################################################