fig2.py 文件源码

def estimate_MSE_vs_alpha(transform, Ms, alphas, K):
    # Without loss of generality
    Z = 1
    tZ = transform(Z)

    # Estimate MSEs by constructing estimators K times
    MSEs = np.empty((len(Ms), len(alphas)))
    MSEs_stdev = np.empty((len(Ms), len(alphas)))
    for Mi, M in enumerate(Ms):
        # Compute means (K x alphas) in a loop, as otherwise
        # this runs out of memory with K = 100,000.
        means = np.empty((K, len(alphas)))
        for ai, alpha in enumerate(alphas):
            Ws = np.power(np.random.exponential(1.0, size=(K, M)), alpha)   # (K, M)
            means[:, ai] = np.mean(Ws, axis=1)
            print_progress('M = %d: done %.0f%%' % (M, 100.0 * (ai+1) / len(alphas)))
        print('')

        g = np.power(gamma(1.0 + alphas), 1.0 / alphas)         # (alphas)
        tZ_hats = transform(g * np.power(means, -1.0/alphas))   # (K, alphas)
        SEs = (tZ_hats - tZ) ** 2                               # (K, alphas)
        MSEs[Mi] = np.mean(SEs, axis=0)                         # (alphas)
        MSEs_stdev[Mi] = np.std(SEs, axis=0) / np.sqrt(K)       # (alphas)

    return MSEs, MSEs_stdev