```python
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
import pandas as pd
from scipy.stats import f, t #use for f and t distribution

from IPython.display import display, Markdown, Math
```


```python
#sympy quicker variable naming
def var(v):
    return sp.symbols(v)

#display math latex
def print_lat(s):
    display(Math(s))
```


```python
class LeastSq:
    def __init__(self, pred, expl, df):
        '''
        Takes pandas dataframe, and makes x, Y based on pred Y,
        and explanatory variables x.
        '''
        self.df = df
        self.pred = pred
        self.expl = expl

        self.x, self.Y = self.gen_x_Y()

        self.dim = [self.x.shape, self.Y.shape]
        self.check_dim()
        
    def gen_x_Y(self):
        '''
        Generates x, Y from pd dataframe
        '''
        res_Y = np.matrix(self.df[self.pred].to_numpy())
        res_x = np.matrix([self.df[i].to_numpy() for i in self.expl])
        
        return res_x, res_Y
        
        
        
    
    #check dimension errors
    def check_dim(self):
        if self.dim[0][1] != self.dim[1][1] or self.dim[1][0] != 1:
            raise ValueError(f"Dimensions are incorrect, x: {self.dim[0]} Y: {self.dim[1]}")
    
    def linear(self):
        '''
        Given x, Y: solve for linear equation
        
        f(x1 ... xn) = A_1*x1 ... An * xn + B
                     = [a1, ... an, 1]x
                     = Ax
                     
                     Ax = f
                     
        for b = f             
        Solved via AtAx = Atb
        '''
        
        A = self.x.T
        A = np.hstack((A, np.ones(self.dim[1]).T))
        
        b = self.Y.T
        
        AtA = np.matmul(A.T, A)
        Atb = np.matmul(A.T, b)
        
        ans = np.linalg.solve(AtA, Atb)
        
        return LQ_out(self.pred, self.expl, df, ans, A)


class LQ_out(LeastSq):
    def __init__(self, pred, expl, df, ans, A):
        '''
        Subclass of LeastSq, handles output results
        '''
        super().__init__(pred, expl, df)
        self.A = A
        
        
        #df = n - k - 1, s.t k = # explanatory variables
        self.df = self.dim[0][1]-self.dim[0][0]-1
        
        #ans vector
        self.ans = ans
        self.ans_dim = ans.shape
        
        #coefficients
        self.coefs = self.get_coefs()
        
        #latex/sympy reporting
        self.func = self.f_sympy()
        self.latex = sp.latex(self.func)
        
        #residual info
        self.resid = self.residuals()
        self.SSE = np.linalg.norm(self.resid.T)**2
        self.RSE = self.get_RSE()
        self.resid_quartiles = self.quartiles(self.resid)
        self.sigma = np.sqrt(self.SSE/self.df)
        
        #SST = SSE + SSM
        self.SST = np.linalg.norm(self.get_SST().T)**2
        self.SSM = self.SST - self.SSE
        
        #R-squared
        self.R_SQ, self.R_SQ_adj = self.get_R_SQ()
        
        #f-statistic
        self.F = self.get_F_stats()
        self.p_value = self.get_p_value()
        
        #t-statistic
        self.SE = self.get_SE()
        self.T = self.get_T_stats()
        self.coef_p_value = self.get_coef_p_value()
        self.signifs = self.signif_codes()
        
    def residuals(self):
        '''
        Given A, x, y, where x is the result of solving least squares, return error
        e = b - Ax
        '''
        
        b = self.Y.T
        
        Ax = np.matmul(self.A, self.ans)
        
        return b-Ax
    
    def dup_rows(self, a, indx, num_dups=1):
        '''
        Duplicates a row in matrix num_dups times
        '''
        return np.insert(a,[indx+1]*num_dups,a[[indx],:],axis=0)
    
    def signif(self, p):
        '''
        Finds significance code of single p-value
        '''
        dic = {"***": (0, 0.001),
               "**": (0.001, 0.01),
               "*": (0.01, 0.05),
               ".": (0.05, 0.1),
               " ": (0.1, 1)}
        
        for key, val in dic.items():
            if val[0] <= p < val[1]:
                return key
            
        return "not a p-value"
        
    
    def signif_codes(self):
        '''
        Generate significance code for all cexplanatory variable coefficients
        
        range based on values in above method dictionary.
        '''
        return [self.signif(i) for i in self.coef_p_value]
        
        
    
    def get_coef_p_value(self):
        '''
        Gets p-value of explanatory variable coefficients s.t:
        
        p_i = 2*(1-T^-1(|t|, df)), i = 1 ... k
        '''
        
        ps = 2*(1-t.cdf(np.abs(self.T), self.df ))
        
        return ps
        
    
    def get_T_stats(self):
        '''
        Gets t statistic s.t:
        
        t_i = beta_i/SE_i, i = 1 ... k
        '''
        
        return self.ans/self.SE.T
    
    def get_SE(self):
        '''
        Gets Standard Error of explanatory variable coefficients s.t:
        
        SE_i = sqrt((AtA)^-1 * S^2), for SEi_i, i = 1 ... k
        '''
        
        AtA_inv = np.linalg.inv(np.matmul(self.A.T, self.A))
        
        AtA_inv_diags = AtA_inv.diagonal()
        
        SE = np.sqrt(AtA_inv_diags*self.sigma**2)
        
        return SE
    
    def get_p_value(self):
        '''
        Gets p-value of f-statistic s.t:
        
        p = 1 - F^-1(f, k, df)
        '''
        x_dim = self.dim[0]
        
        p = 1 - f.cdf(self.F, x_dim[0], self.df)
        
        return float("%.4g" % p)
        
    
    def get_F_stats(self):
        '''
        Gets F-statistic s.t:
        
        f = (SSM/k)/(SSE/df)
        '''
        x_dim = self.dim[0]
        F = (self.SSM/x_dim[0])/(self.SSE/(self.df))
        
        return round(F, 2)
    
    def get_R_SQ(self):
        '''
        Gets both R-squared, and Adjusted R-squared s.t:
        
        R-squared = SSM/SST
        
        Adj. R-squared: (n-1)/df * R-squared - k/df
        '''
        rsq = self.SSM/self.SST
        
        x_dim = self.dim[0]
        rsq_adj = ((x_dim[1]-1)/(self.df))*rsq - (x_dim[0])/(self.df)
        
        return np.round([rsq, rsq_adj], 4)
    
    def get_SST(self):
        '''
        Gets SST, the square difference between each predictor variable
        and the mean of the predictor s.t:
        
        SST = sum (y_i - ybar)^2
        '''
        b = self.Y.T
        
        y_bar = np.array([[np.mean(self.Y)]]*self.dim[0][1])
        
        return b-y_bar
    
    def quartiles(self, dat):
        '''
        Finds quantiles, and min/max based on sorted dat, (meant for residuals)
        '''
        sorted_dat = np.sort(dat.A1)
        
        Qs_tags = ["Min", "1Q", "Median", "3Q", "Max"]
        
        Qs = {i : round(np.quantile(sorted_dat, indx*.25), 3) for indx, i in enumerate(Qs_tags)}

        return Qs
        
    def get_RSE(self):
        '''
        Gets Residual standard error, s.t:
        
        RSE = sqrt(SSE/df)
        '''
        
        RSE = np.sqrt(self.SSE/self.df)
        
        return round(RSE, 1)
        
    def get_coefs(self):
        '''
        Gets coefficients as a dictionary based on explantory variables,
        and intercept
        
        s.t linear equation is: tags["Intercept"] + sum (tags[expl] * x_i)
        '''
        tags = self.expl + ["Intercept"]
        return {i : round(float(self.ans[indx]), 3) for indx, i in enumerate(tags)}
    
    def f_sympy(self):
        '''
        Generate sympy function based on result of regression model.
        
        variable names based on chosen explanatory variables, in paranthesis.
        '''
        eq = 0
        
        for i in self.expl:
            eq += self.coefs[i] * var(f"({i})")
            
        eq += self.coefs["Intercept"]
        
        return eq
    
    def summary(self):
        '''
        Does summary function using computed values, based on the
        R function summary() for regression models.
        '''
        
        #Do Call
        sum_expls = " + ".join(x for x in self.expl)
        print(f"Call:\nlm(formula = {self.pred} ~ {sum_expls}, data = df)")
        
        #Do Residuals
        resid_tags = "\tMin \t 1Q\t Median\t 3Q \t Max"
        resid_qs = "\t".join(str(x) for x in list(self.resid_quartiles.values()))
        print("Residuals:")
        print(resid_tags.expandtabs(4))
        print(resid_qs, "\n")
        
        #Do Coefficients
        expl_tags = ["(Intercept)"] + self.expl
        #print(self.signifs)
        
        dats = [self.ans, self.SE.T, self.T, self.coef_p_value, ["{:<4}".format(i) for i in self.signifs]]
        dats = [np.asarray(i).reshape(-1) for i in dats]
    
    
        df = pd.DataFrame()
        df["Estimate"] = dats[0]
        df["Std. Error"] = dats[1]
        df["t value"] = dats[2]
        df["Pr(>|t|)"] = dats[3]
        df[""] = dats[4]
        

        df = df.apply(np.roll, shift=1)
        
        df.index = expl_tags
        
        print("Coefficients:")
        print(df)
        
        
        #Do Signif
        print("---")
        print("Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1\n")
        
        #Do Resid SE
        print("Residual standard error:", self.RSE, "on", self.df,  "degrees of freedom")
        
        #Do Multiple R-Sq
        print(f"Multiple R-squared: {self.R_SQ}, Adjusted R-squared: {self.R_SQ_adj}")
        
        #Do F-stat
        print(f"F-statistic: {self.F}, on {self.dim[0][0]} and {self.df} DF, p-value: {self.p_value}")

            

        
        
        
```


```python
#Example (requires GPA.csv to work)
df = pd.read_csv("data/GPA.csv")
df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>obs</th>
      <th>GPA</th>
      <th>HSM</th>
      <th>HSS</th>
      <th>HSE</th>
      <th>SATM</th>
      <th>SATCR</th>
      <th>SATW</th>
      <th>sex</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1</td>
      <td>3.84</td>
      <td>10</td>
      <td>10</td>
      <td>10</td>
      <td>630</td>
      <td>570</td>
      <td>590</td>
      <td>2</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2</td>
      <td>3.97</td>
      <td>10</td>
      <td>10</td>
      <td>10</td>
      <td>750</td>
      <td>700</td>
      <td>630</td>
      <td>1</td>
    </tr>
    <tr>
      <th>2</th>
      <td>3</td>
      <td>3.49</td>
      <td>8</td>
      <td>10</td>
      <td>9</td>
      <td>570</td>
      <td>510</td>
      <td>490</td>
      <td>2</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4</td>
      <td>1.95</td>
      <td>6</td>
      <td>4</td>
      <td>8</td>
      <td>640</td>
      <td>600</td>
      <td>610</td>
      <td>1</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5</td>
      <td>2.59</td>
      <td>8</td>
      <td>10</td>
      <td>9</td>
      <td>510</td>
      <td>490</td>
      <td>490</td>
      <td>2</td>
    </tr>
  </tbody>
</table>
</div>




```python
model = LeastSq("GPA", ["HSM", "HSS", "HSE"], df)
```


```python
res = model.linear()
```


```python
res.summary()
```

    Call:
    lm(formula = GPA ~ HSM + HSS + HSE, data = df)
    Residuals:
        Min      1Q  Median  3Q      Max
    -2.203	-0.393	0.142	0.545	1.496 
    
    Coefficients:
                 Estimate  Std. Error   t value  Pr(>|t|)      
    (Intercept)  0.069304    0.453656  0.152767  0.878793      
    HSM          0.123246    0.054855  2.246774  0.026155  *   
    HSS          0.136137    0.069951  1.946177  0.053554  .   
    HSE          0.058478    0.065417  0.893916  0.372838      
    ---
    Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
    
    Residual standard error: 0.7 on 146 degrees of freedom
    Multiple R-squared: 0.2277, Adjusted R-squared: 0.2119
    F-statistic: 14.35, on 3 and 146 DF, p-value: 3.036e-08



```python

```

2c2a133ec58ecafdfd2d42662ae04214b1f314fc
ipynb
Jupyter Notebook
Least Squares.ipynb
AlephEleven/Python-lm-function
e5f0bb42e9357a29d04f8a6a785538ef7c5a398a
MIT
Least Squares.ipynb
AlephEleven/Python-lm-function
e5f0bb42e9357a29d04f8a6a785538ef7c5a398a
MIT
Least Squares.ipynb
AlephEleven/Python-lm-function
e5f0bb42e9357a29d04f8a6a785538ef7c5a398a
MIT
Qwen/Qwen-72B
1. YES
2. YES