OLSの検定・区間推定

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf

beta = [10, 5]

def gen_data(n = 300, random_state=0):
    np.random.seed(random_state)
    x = np.random.uniform(size=n)
    y = beta[0] + x * beta[1] + np.random.normal(size=n)
    return x, y

estimates = []
for seed in range(1000):
    x, y = gen_data(n = 100, random_state=seed)
    df = pd.DataFrame(dict(x=x, y=y))
    beta1 = sum((x - x.mean()) * (y - y.mean())) / sum((x - x.mean())**2)
    estimates.append(beta1)
estimates = np.array(estimates)


fig, ax = plt.subplots()
ax.hist(estimates, bins=20)
ax.axvline(beta[1], color="gray", label=r"True $\beta_1$")
ax.set(title=r"Distribution of $\hat{\beta}_1$", xlabel=r"$\hat{\beta}_1$")
ax.legend()

from myst_nb import glue
glue("dist_beta_fig", fig, display=False)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import norm


beta = [10, 5]


def gen_data(n = 300, random_state=0):
    np.random.seed(random_state)
    x = np.random.uniform(size=n)
    y = beta[0] + x * beta[1] + np.random.normal(size=n)
    return x, y


def se_beta(x, residuals):
    n = len(x)
    sqd = (x - np.mean(x)) ** 2
    df = 2
    nu = (1 / (n - df)) * np.sum( sqd * (residuals ** 2) )
    de = ( (1 / n) * np.sum(sqd) ) ** 2
    var = (1 / n) * (nu) / (de)
    return np.sqrt(var)


fig, ax = plt.subplots()

estimates = []
se_list = []
for seed in range(100):
    x, y = gen_data(n = 300, random_state=seed)
    df = pd.DataFrame(dict(x=x, y=y))
    model = smf.ols('y ~ x', data=df).fit()

    beta1 = model.params["x"]
    se_beta1 = se_beta(x, residuals=model.resid)

    x_plot = np.linspace(4, 6, 50)
    y_plot = norm.pdf(x=x_plot, loc=beta1, scale=se_beta1)
    ax.plot(x_plot, y_plot, alpha=.5)

    estimates.append(beta1)
    se_list.append(se_beta1)
estimates = np.array(estimates)


ax.hist(estimates, bins=10)
ax.axvline(beta[1], color="gray", label=r"True $\beta_1$")
ax.set(title=r"Distribution of $\hat{\beta}_1$", xlabel=r"$\hat{\beta}_1$")
ax.legend()

from myst_nb import glue
glue("dist_beta_fig", fig, display=False)

# import numpy as np
# import pandas as pd

# n = 300
# beta = [10, 5]
# np.random.seed(0)
# x = np.random.uniform(size=n)
# df = pd.DataFrame({
#     "x": x,
#     "y": beta[0] + x * beta[1] + np.random.normal(size=n)
# })
# # df.plot.scatter(x="x", y="y")

# import statsmodels.api as sm
# import statsmodels.formula.api as smf
# from stargazer.stargazer import Stargazer

# models = [
#     smf.ols('y ~ x', data=df).fit(),
# ]
# Stargazer(models)