<jupyter_start><jupyter_text>Simulations for Problem 1<jupyter_code>import collections
import itertools

import matplotlib.pyplot as plt
import multiprocessing
import numpy as np
import pandas as pd
from scipy import linalg
from scipy import stats

X_VALUES = [0, 1, 3]
MU = np.array([0.6, 3.6, np.nan, 6.8])


def simulate_data(n, random_seed=None):
    np.random.seed(random_seed)
    dist = stats.rv_discrete(values=list(zip(
        *[(x, 1/len(X_VALUES)) for x in X_VALUES])))
    x = dist.rvs(size=(n,))
    while len(np.unique(x)) == 1:
        x = dist.rvs(size=(n,))
    mu = MU[x]
    return x, stats.norm.rvs(loc=mu, scale=1)

def solve_beta_hat(x, y):
    X = np.column_stack((np.ones_like(x), x))
    return linalg.cho_solve(linalg.cho_factor(X.T.dot(X)), np.eye(X.shape[1])).dot(X.T).dot(y)
x, y = simulate_data(20, random_seed=2019)
beta_hat = solve_beta_hat(x, y)
beta_hat
fig = plt.figure(figsize=(6,4))
ax = fig.gca()
ax.grid(True)
ax.plot(x, y, 'o')
xlim_left, xlim_right = ax.get_xlim()
ax.plot([xlim_left, xlim_right],
        np.array([xlim_left, xlim_right])*beta_hat[1] + beta_hat[0], '--k')
ax.set_xlim(xlim_left, xlim_right);
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_title('Ordinary Least Squares Regression')
fig.tight_layout()
fig.savefig('p1_ols.pdf', bbox_inches='tight')
def _estimate_beta1(n):
    x, y = simulate_data(n)    
    return solve_beta_hat(x, y)[1]

def estimate_beta1(num_estimates, n):    
    pool = multiprocessing.Pool(4)
    estimates = pool.map(_estimate_beta1, [n]*num_estimates)
    pool.close()
    return np.array(estimates)

def estimate_bias(num_trials):
    n_list = list(range(2, 31, 2)) + list(range(30, 101, 5))
    mean = []
    standard_errors = []
    for n in n_list:
        estimates = estimate_beta1(num_trials, n)
        mean.append(np.mean(estimates))
        standard_errors.append(np.sqrt(np.var(estimates, ddof=1)))
    return pd.DataFrame(collections.OrderedDict([
        ('$n$', n_list),
        ('$\\hat{\\beta}_1$', mean),
        ('$\\hat{\\sigma}$', standard_errors),
    ]))

#bias_simulations = estimate_bias(1000000)
#bias_simulations.to_pickle('p1_simulations.pkl')
bias_simulations = pd.read_pickle('p1_simulations.pkl')
bias_simulations
def partition(n, k = None, m = None):
    k = n if k is None else k
    m = n if m is None else m
    if m is None or m >= n:
        yield [n]
        
    for f in range(n-1 if m >= n else m, 0, -1):
        if f*(k - 1) < n - f:
            break
        for p in partition(n-f, k - 1, f):
            yield [f] + p
            

def multiset_permuation(xs):
    current_permutation = len(xs)*[None]
    counts = collections.Counter(xs)
    
    def yield_permuation(j):
        if j == len(xs):
            yield tuple(current_permutation)
            return
            
        for x in counts:
            if counts[x] <= 0:
                continue                
            counts[x] -= 1
            current_permutation[j] = x
            for permuation in yield_permuation(j + 1):
                yield permuation
            counts[x] += 1
    
    return yield_permuation(0)
            
def enumerate_multinomial_support(n, k):
    for p in partition(n, k):
        unordered_draw = [0]*(k - len(p)) + p
        for x in multiset_permuation(unordered_draw):
            yield x
def solve_expected_beta(n):
    support = enumerate_multinomial_support(n, 3)
    weighted_sum = np.zeros(2)
    total_prob = 0
    for sample in support:
        if np.isin(n, sample):
            continue        
        x = np.hstack([count*[X_VALUES[i]] for i, count in enumerate(sample)])
        y = MU[x.astype(np.int)]
        expected_beta = solve_beta_hat(x, y)
        prob = stats.multinomial.pmf(sample, n, p = np.ones(len(X_VALUES))/len(X_VALUES))
        weighted_sum += expected_beta*prob
        total_prob += prob
    return weighted_sum/total_prob

expected_beta = np.array([solve_expected_beta(n) for n in bias_simulations['$n$']])
expected_beta
bias_simulations['$\\mathbb{E}[\\hat{\\beta}_1]$'] = expected_beta[:,1]
bias_simulations
with open('p1_simulation_results.tex', 'w') as f:
    f.write(bias_simulations.to_latex(escape=False, index=False))
fig = plt.figure(figsize=(6, 4))
ax = fig.gca()
ax.grid(True)

ax.plot(bias_simulations['$n$'], bias_simulations['$\\hat{\\beta}_1$'],
        '.-', label='$\\hat{\\beta}_1$', c=plt.cm.Set1(1))
ax.plot(bias_simulations['$n$'], bias_simulations['$\\hat{\\sigma}$'],
        '.-', label='$\\hat{\\sigma}$', c=plt.cm.Set1(2))
ax.plot(bias_simulations['$n$'], bias_simulations['$\\hat{\\beta}_1$'] - 2,
        '.-', label='Bias($\\hat{\\beta}_1$)', c=plt.cm.Set1(3))
ax.plot(bias_simulations['$n$'], bias_simulations['$\\hat{\\sigma}$']**2,
        '.-', label='$\\hat{\\sigma}^2$', c=plt.cm.Set1(4))
ax.hlines(y=2, xmin=ax.get_xlim()[0], xmax=ax.get_xlim()[1],
          linestyles='--', label='$\\beta_1$')
ax.legend()
ax.set_xlabel('$n$')
ax.set_title('Simulation Results')
fig.tight_layout()
fig.savefig('p1_simulation_results.pdf', bbox_inches='tight')
fig = plt.figure(figsize=(6, 4))
ax = fig.gca()
ax.grid(True)

ax.plot(bias_simulations['$n$'], bias_simulations['$\\hat{\\beta}_1$'] - 2,
        '.-', label='Simulated Bias($\\hat{\\beta}_1$)', c=plt.cm.Set1(0), alpha=0.4)
ax.plot(bias_simulations['$n$'], bias_simulations['$\mathbb{E}[\\hat{\\beta}_1]$'] - 2,
        '.-', label='Bias($\\hat{\\beta}_1$)', c=plt.cm.Set1(1), alpha=0.4)
ax.legend()
ax.set_title('Bias Comparison')
ax.set_xlabel('$n$')
ax.set_ylabel('Bias')
ax.set_yscale('symlog', linthreshy=1e-5)
fig.tight_layout()
fig.savefig('p1_bias_comparison.pdf', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17952783650458068, 'fraction_numerical': 0.01937984496124031, 'mean_word_length': 3.2877643504531724, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '180964', 'license': 'MIT', 'n_tokens_mistral': 2290, 'n_tokens_neox': 2217, 'n_words': 396, 'path': 'hw1/p1_simulation.ipynb', 'repo_name': 'ppham27/stat571', 'size': '5676'}