import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import statsmodels.api as sm
import scipy.stats as stats

n = 100
p = 3
np.random.seed(201681)
X1 = np.random.chisquare(5, n)
X2 = np.random.chisquare(3, n)
X3 = np.random.uniform(3, 9, n)
E = np.random.normal(5, 7, n)
Y = 7 + 3*X1 + 7*X2 + 5*X3 + E

data = np.column_stack((Y, X1, X2, X3))
DT = pd.DataFrame(data, columns=["Y", "X1", "X2", "X3"])
DT.head()
DT.tail()

# Bootstrap Sample
coef = np.zeros((n, p + 1))

for i in range(n):
    ID = np.random.choice(DT.index, n, replace = True)
    sample = DT.iloc[ID, :]
    X_i = sample.drop(["Y"], axis = 1)
    X_i = sm.add_constant(X_i) 
    Y_i = sample["Y"]
    model = sm.OLS(Y_i, X_i)
    result = model.fit()
    coef[i, 0] = result.params[0]
    coef[i, 1] = result.params[1]
    coef[i, 2] = result.params[2]
    coef[i, 3] = result.params[3]

beta = pd.DataFrame(coef, columns=["B0", "B1", "B2", "B3"])

# Plotting
plt.figure(figsize = (12, 12))
#1
plt.subplot(4, 2, 1)
sns.histplot(beta["B0"], bins = 20, kde = True, color = "blue" )
plt.xlabel("Value of B0")
plt.ylabel("Frequency")
plt.title("Histogram and Density Plot for B0")

plt.subplot(4, 2, 2)
stats.probplot(beta["B0"], dist="norm", plot = plt)
plt.title("QQ Plot for B0")

#2
plt.subplot(4, 2, 3)
sns.histplot(beta["B1"], bins = 20, kde = True, color = "red" )
plt.xlabel("Value of B1")
plt.ylabel("Frequency")
plt.title("Histogram and Density Plot for B1")

plt.subplot(4, 2, 4)
stats.probplot(beta["B1"], dist="norm", plot = plt)
plt.title("QQ Plot for B1")

#3
plt.subplot(4, 2, 5)
sns.histplot(beta["B2"], bins = 20, kde = True, color = "green" )
plt.xlabel("Value of B2")
plt.ylabel("Frequency")
plt.title("Histogram and Density Plot for B2")

plt.subplot(4, 2, 6)
stats.probplot(beta["B2"], dist="norm", plot = plt)
plt.title("QQ Plot for B2")

#4
plt.subplot(4, 2, 7)
sns.histplot(beta["B3"], bins = 20, kde = True, color = "blueviolet" )
plt.xlabel("Value of B3")
plt.ylabel("Frequency")
plt.title("Histogram and Density Plot for B3")

plt.subplot(4, 2, 8)
stats.probplot(beta["B3"], dist="norm", plot = plt)
plt.title("QQ Plot for B3")

plt.tight_layout()
plt.savefig("D:/Histograms and Q-Q Plots.png")
