<jupyter_start><jupyter_text>Utils<jupyter_code># Gauss 
def gauss(x,mu,s):
    #:params x: 1-D array data 
    #:prams mu,s: \mu and \sigma 
    #:return: \phi(x) 
    d = x.shape[0]
    m = mu.shape[0]
    phi = np.exp(-(x.reshape(-1,1) - mu)/(2*s**2))
    phi = np.concatenate(([1],phi.ravel()))
    return phi

# Function
def f(x):
    return 3*np.sin(x)

# Make Toy Data 
def make_toy_data(N,lower = 0,upper = 2*np.pi,std = 1):
    X = np.random.rand(N)*(upper-lower) + lower
    y = f(X) + np.random.randn(N)*std 
    return X.reshape(-1,1),y.reshape(-1,1)


# Color maps 
cmaps = [[0.122, 0.467, 0.706],"orange","green"]

# Plot Prediction 
def plot_prediction(X_tr,y_tr,regressor,title,lower = 0,upper = 2*np.pi):
    X = np.linspace(lower,upper,100).reshape(-1,1)
    y_pred = regressor.predict(X) 
    y_true = f(X)
    
    rmse = np.mean((y_true-y_pred)**2)**0.5
    print(f"RMSE : {rmse}")
    
    fig,ax = plt.subplots(1,1,figsize = (10,7))
    ax.plot(X,y_pred,label="Predict",color=cmaps[0])
    ax.plot(X,y_true,label="Ground Truth",color=cmaps[1])
    ax.scatter(X_tr,y_tr,label="Training Data",color=cmaps[2])
    ax.set_title(title)
    
    plt.legend()
    plt.show()
    

# Plot Prediction with Std 
def plot_prediction_with_std(X_tr,y_tr,bayes_regressor,title,lower = 0,upper = 2*np.pi):
    X = np.linspace(lower,upper,100).reshape(-1,1)
    y_pred,y_std = bayes_regressor.predict(X,return_std=True)
    y_true = f(X)
    
    rmse = np.mean((y_true-y_pred)**2)**0.5
    print(f"RMSE : {rmse}")
    
    fig,ax = plt.subplots(1,1,figsize = (10,7))
    ax.plot(X,y_pred,label="Predict",color=cmaps[0])
    
    y_pred_upper = y_pred + y_std
    y_pred_lower = y_pred - y_std 
    ax.fill_between(X.ravel(),y_pred_lower.ravel(),y_pred_upper.ravel(),alpha=0.3,color=cmaps[0])
    
    ax.plot(X,y_true,label="Ground Truth",color=cmaps[1])
    ax.scatter(X_tr,y_tr,label="Training Data",color=cmaps[2])
    ax.set_title(title)
    
    plt.legend()
    plt.show()<jupyter_output><empty_output><jupyter_text>Maximum Likelihood where $\Phi$ is a design matrix, $T$ is target matrix      $$W_{ML} = (\Phi^T\Phi)^{-1}\Phi^T T$$      $$\frac{1}{\beta_{ML}} = \frac{1}{NK} \sum_{n = 1}^N ||t_n - W^T\phi(x_n)||^2$$<jupyter_code>class LinearRegression():
    def __init__(self,mu = None,s = None):
        self.weight = None
        self.beta = None 
        self.mu = mu
        self.s = s
        self.phi = lambda x:gauss(x,self.mu,self.s)
        
    def fit(self,X,y):
        #:params X: 2-D array (N_samples,N_dims)
        #:params y: 2-D array (N_samples,N_targets)  
        N = X.shape[0] 
        K = y.shape[1] 
        design_mat = np.vstack([self.phi(x) for x in X])
        self.weight = np.linalg.inv(design_mat.T@design_mat)@design_mat.T@y  
        tmp = y - self.weight.T@design_mat.T 
        self.beta = N*K/np.sum(tmp**2) 
        
    def predict(self,X):
        #:params X: 2-D array (N_samples,N_dims) N_dims = len(mu) = len(s) 
        design_mat = np.vstack([self.phi(x) for x in X])
        return np.dot(design_mat,self.weight)
X_tr,y_tr = make_toy_data(100,std = 0.75)
mu = np.random.rand(30)*6 - 3
s = np.random.rand(30)*5
lr = LinearRegression(mu = mu,s = s)
lr.fit(X_tr,y_tr)
plot_prediction(X_tr,y_tr,lr,"LinearRegression")<jupyter_output>RMSE : 0.9931091929541018<jupyter_text>Ridge    almost same as LinearRegression but it has a regularization term.<jupyter_code>class Ridge():
    def __init__(self,lamda=1e-2,mu = None,s = None):
        self.weight = None
        self.lamda = lamda 
        self.mu = mu
        self.s = s
        self.phi = lambda x:gauss(x,self.mu,self.s)
        
    def fit(self,X,y):
        #:params X: 2-D array (N_samples,N_dims)
        #:params y: 2-D array (N_samples,N_targets) 
        M = X.shape[1]*self.mu.shape[0] + 1
        design_mat = np.vstack([self.phi(x) for x in X])
        self.weight = np.linalg.inv(self.lamda*np.eye(M) + design_mat.T@design_mat)@design_mat.T@y 

    def predict(self,X):
        #:params X: 2-D array (N_samples,N_dims) N_dims = len(mu) = len(s) 
        design_mat = np.vstack([self.phi(x) for x in X])
        return np.dot(design_mat,self.weight)
ridge = Ridge(lamda = 1e-8,mu = mu,s = s)  
ridge.fit(X_tr,y_tr)

plot_prediction(X_tr,y_tr,ridge,"Ridge")<jupyter_output>RMSE : 0.23627845905908115<jupyter_text>Bayesian Linear Regression  Prior Distributin for Weight $$p(w|\alpha) = \mathcal{N}(w|0,\alpha^{-1}I)$$  Posterior Distribution for Weights $$p(w|t) = \mathcal{N}(w|m_N,S_N)$$ where $$m_N = \beta S_N \Phi^\top t$$ $$S_N^{-1} = \alpha I + \beta \Phi^\top \Phi$$  Predictive Distribution  $$p(t|x) = \mathcal{N}(t|m_N^\top\phi(x),\sigma_N^2(x))$$ where $$\sigma_N^2(x) = \frac{1}{\beta} + \phi(x)^\top S_N\phi(x)$$  Evidence  $$\ln{p(t|\alpha,\beta)} = \frac{M}{2}\ln{\alpha} + \frac{N}{2}\ln{\beta} - E(m_N) + \frac{1}{2}\ln{|S_N|} - \frac{N}{2}\ln{2\pi}$$ where $$E(m_N) = \frac{\beta}{2}||t - \Phi m_N||^2 + \frac{\alpha}{2}m_N^\top m_N$$  Maximizing the evidence function Let $\lambda_i$ be eigenvalue of $\beta\Phi^\top\Phi$    $$\gamma = \sum_{i = 1}^M \frac{\lambda_i}{\alpha + \lambda_i}$$      $$\alpha = \frac{\gamma}{m_N^\top m_N}$$      $$\frac{1}{\beta} = \frac{1}{N - \gamma}||t - \Phi m_N||^2$$<jupyter_code>class BayesianLinearRegression():
    def __init__(self,alpha = 1e-1,beta = 1e-1,mu = None,s = None):
        self.weight = None
        self.S = None
        self.M = None
        self.N = 0
        self.alpha = alpha 
        self.beta = beta 
        self.mu = mu
        self.s = s
        self.phi = lambda x:gauss(x,mu,s)
        
    def fit(self,X,y,optimize_evidence = False,n_iters = 20,threshold = 1e-3):
        #:params X: 2-D array (N_samples,N_dims)
        #:params y: 1-D array (N_samples)  
        #:params optimze_evidence: if alpha and beta is optimized or not 
        
        self.N = X.shape[0]
        self.M = X.shape[1]*self.mu.shape[0] + 1
        
        if optimize_evidence:
            self.optimize_evidence_(X,y,n_iters,threshold) 
            
        design_mat = np.vstack([self.phi(x) for x in X])
        self.S = np.linalg.inv(self.alpha*np.eye(self.M) + self.beta*design_mat.T@design_mat) 
        self.weight = self.beta*self.S@design_mat.T@y 
    
    def partial_fit(self,X,y):
        # Before this method is called, fit() should be called 
        
        self.N += X.shape[0] 
        design_mat = np.vstack([self.phi(x) for x in X])
        S_old_inv = np.linalg.inv(self.S) 
        
        self.S = np.linalg.inv(S_old_inv + self.beta*design_mat.T@design_mat)
        self.weight = self.S@(S_old_inv@self.weight + self.beta*design_mat.T@y)
    
    def calc_evidence_(self,tmp):
        E = self.beta/2*tmp + self.alpha/2*np.dot(self.weight,self.weight) 
        evidence = self.M*np.log(self.alpha)/2 + self.N*np.log(self.beta)/2 - E + np.linalg.det(self.S) - self.N*np.log(2*np.pi)/2 
        return evidence 
    
    def optimize_evidence_(self,X,y,n_iters,threshold):
        #:params n_iters: Number of times to optimize alpha and beta 
        #:params threshold: If the difference of evidence is lower than this, 
        
        design_mat = np.vstack([self.phi(x) for x in X])
        C = design_mat.T@design_mat
        org_lambdas,_ = np.linalg.eig(C)
        with warnings.catch_warnings(): # Ignore Warnings
            warnings.simplefilter('ignore')
            org_lambdas = org_lambdas.astype(np.float64) 
        before_evidence = -10**10 
        
        for _ in range(n_iters):
            self.S = np.linalg.inv(self.alpha*np.eye(self.M) + self.beta*C) 
            self.weight = self.beta*self.S@design_mat.T@y 
            
            lambdas = self.beta*org_lambdas
            gamma = np.sum(lambdas/(lambdas + self.alpha))
            self.alpha = gamma/np.dot(self.weight,self.weight)
            tmp = y - design_mat@self.weight
            tmp = np.dot(tmp,tmp)
            self.beta = (self.N - gamma)/tmp 
            evidence = self.calc_evidence_(tmp)
            
            if np.abs(before_evidence-evidence) < threshold:
                break
            before_evidence = evidence
        
    def predict(self,X,return_std = False):
        design_mat = np.vstack([self.phi(x) for x in X])
        pred = np.dot(design_mat,self.weight).ravel()
        if return_std:
            std = np.sqrt(1/self.beta + np.diag(design_mat@self.S@design_mat.T))
            return pred,std
        else:
            return pred
blr = BayesianLinearRegression(alpha = 1e-7,beta = 1.3,mu = mu,s = s)
blr.fit(X_tr,y_tr.ravel())

plot_prediction_with_std(X_tr,y_tr,blr,"BayesianLinearRegression")
# partial fit

X_new,y_new = make_toy_data(100,lower = 0,upper = np.pi,std = 0.5)

blr.partial_fit(X_new,y_new.ravel())
plot_prediction_with_std(np.concatenate([X_tr,X_new]),
                         np.concatenate([y_tr,y_new]),
                         blr,"BayesianLinearRegression",lower = 0,upper = 2*np.pi)
# With Evidence Optimization 

blr = BayesianLinearRegression(alpha = 10,beta = 10,mu = mu,s = s)
blr.fit(X_tr,y_tr.ravel(),optimize_evidence = True,n_iters = 100,threshold = 1e-8)
print(f"Optimized alpha : {blr.alpha}")
print(f"Optimized beta : {blr.beta}")

plot_prediction_with_std(X_tr,y_tr,blr,"BayesianLinearRegression with EvidenceOptimization")<jupyter_output>Optimized alpha : 4.15731245173827e-06
Optimized beta : 2.0158980963232715
RMSE : 3.058773617893491
{'fraction_non_alphanumeric': 0.16300307953700754, 'fraction_numerical': 0.0270786874800892, 'mean_word_length': 2.779293739967897, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '458835', 'license': 'MIT', 'n_tokens_mistral': 3949, 'n_tokens_neox': 3704, 'n_words': 617, 'path': 'notebook/chapter03_linear_regression.ipynb', 'repo_name': 'hedwig100/PRML', 'size': '9417'}