<jupyter_start><jupyter_text>Define a random polynomial10th order.<jupyter_code>torch.manual_seed(42)
random_coeffs = torch.randn(10)
random_poly = np.polynomial.polynomial.Polynomial(random_coeffs)
random_poly
#plt.plot(*random_poly.linspace(100, [-5,5]))<jupyter_output><empty_output><jupyter_text>Compute the actual expectation under $p^*(x)$ for bias calculations. Use quadrature.<jupyter_code>random_poly_expectation, error_est = quad(lambda x: log_p(x*torch.ones(1)).exp() * random_poly(x),
                              -15,
                               15,
                               limit=1000)
random_poly_expectation, error_est<jupyter_output>/tmp/ipykernel_19849/3063319930.py:1: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  random_poly_expectation, error_est = quad(lambda x: log_p(x*torch.ones(1)).exp() * random_poly(x),<jupyter_text>Bias, variance as a function of $\lambda$ Pick a few values of $\lambda$.<jupyter_code># NEED TO GET THESE RANGES BETTER
r=15
configs = [{'lambda': 1.0, 'log_sig_range': (-16, -2), 'mu_range': (-r, r)},]
#            {'lambda': 2.0, 'log_sig_range': (-3, 1.5), 'mu_range': (-r, r)},
#            {'lambda': 8.0, 'log_sig_range': (-2, 1.5), 'mu_range': (-r, r)},
#            {'lambda': 64.0, 'log_sig_range': (-1, 1.5), 'mu_range': (-r, r)}]
for l in np.logspace(0,2,20,)[1:]:
    configs.append({'lambda': l, 'log_sig_range': (-5, 1.5), 'mu_range': (-r, r)})<jupyter_output><empty_output><jupyter_text>Now calculate bias and variance<jupyter_code>def sample_theta(psi, n, mm, gg, and_idx=False):
    idx = torch.multinomial(psi.flatten(), num_samples=n, replacement=True)
    if not and_idx: return torch.stack([mm.flatten()[idx], gg.flatten()[idx]], dim=1)
    return torch.stack([mm.flatten()[idx], gg.flatten()[idx]], dim=1), idx
    

def log2prob(logp):
    p = torch.exp(logp - logp.max())
    return p / p.sum()

def log_q(x, theta):
    sig2 = torch.exp(2*theta[1])
    return -0.5*(x.view(-1,1)-theta[0])**2/sig2 - log(2*pi)/2 - theta[1]
def MC_integrate(conf, Neff, fn, plot=False):
    """## Self-normalizing importance sampling integration
    # Propose samples of theta from a big Gaussian or uniform
    # For each theta we calculate log psi with Stam's
    # Then weight the contribution of the integral over that kernel by the proposal dist
    # iterate until the effective sample size reaches Neff
    """
    # todo don't preallocate samples; do online. legacy from before ESS while-loop architecture
    if conf['lambda']>1.1:
        #from a normal dist?
        proposal_p = torch.distributions.multivariate_normal.MultivariateNormal(torch.zeros(2), 2*torch.eye(2))
        samples = proposal_p.sample([100*Neff])
    else:
        # or from a uniform dist in a box supplied by the ranges?
        proposal_p=None
        samples = torch.rand(100*Neff,2)
        samples[:,0] = samples[:,0]*(conf['mu_range'][1] - conf['mu_range'][0]) + conf['mu_range'][0]
        samples[:,1] = samples[:,1]*(conf['log_sig_range'][1] - conf['log_sig_range'][0]) + conf['log_sig_range'][0]

    q = MVNIso(d=1)
    integral = 0
    curve = torch.zeros_like(x_grid)
    weights = []
    weights2 = []
    sum_weights = 0
    sum_weights2 = 0
    ESS = 0
    i = 0
    while ESS <= Neff and i < 100*Neff:
        sample = samples[i]
        kernel_exp = q.set_theta(sample).quadrature_ev(fn)
        log_psi = get_log_psi(q, conf['lambda'])
        log_prop = 0 if proposal_p is None else proposal_p.log_prob(sample)
        discount = (log_psi - log_prop).exp()
        weights.append(discount)
        weights2.append(discount**2)
        
        if torch.isnan(discount) or torch.isnan(kernel_exp): 
            continue
            
        integral = integral + discount * kernel_exp
        sum_weights += discount
        sum_weights2 += weights2[-1]
        ESS = sum_weights**2 / sum_weights2
        i+=1
        
        if plot: curve += discount * q.to_torch_mvn().log_prob(x_grid[:,None]).exp()

    
    if plot:
        plt.plot(x_grid, curve/sum_weights)
        plt.plot(x_grid, log_p(x_grid).exp(),'r')
        plt.title("EV = {}".format(integral/sum_weights))
        plt.show()
    print("N {} Effective N: {}".format(i,ESS.item()))
    return integral / sum_weights
exps_MC = [MC_integrate(conf, 1000, random_poly, True) for conf in configs]
plt.loglog(np.logspace(0,1,10,),np.abs(np.array(exps_MC)),'o-')
plt.hlines(random_poly_expectation, 1,100,color='r')
plt.xlabel("$\lambda$");
plt.ylabel("EV")
plt.show()<jupyter_output><empty_output><jupyter_text>See the convergence towards a single value with increasing N<jupyter_code>exps_MC_conv = [MC_integrate(configs[2], int(n), random_poly) for n in np.logspace(2,4,10)]
plt.plot(exps_MC_conv)<jupyter_output><empty_output><jupyter_text>Now check out the variance as a function of $\lambda$<jupyter_code>def get_mean_sq(conf, fn, n_runs = 50, N=100):
    sq_diffs = []
    for i in range(n_runs):
        estimate = MC_integrate(conf, N, fn)
        sq_diffs.append((estimate)**2)
    return np.mean(sq_diffs)
variances = np.array([get_mean_sq(c,random_poly) for c in configs]) - np.array(exps_MC)**2
plt.semilogx(np.logspace(0,2,20,),np.abs(np.array(exps_MC)-random_poly_expectation))
# plt.plot(variances)
plt.xlabel("$\lambda$");
plt.ylabel("Bias")
plt.show()
plt.semilogx(np.logspace(0,2,20,),variances)
# plt.plot(variances)
plt.xlabel("$\lambda$");
plt.ylabel("Variance")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15374955148905634, 'fraction_numerical': 0.02834589163975601, 'mean_word_length': 3.7125950972104818, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '592032', 'license': 'MIT', 'n_tokens_mistral': 2101, 'n_tokens_neox': 2015, 'n_words': 417, 'path': 'notebooks/bias_variance.ipynb', 'repo_name': 'wrongu/sampling-variational', 'size': '5574'}