<jupyter_start><jupyter_text>Least squares regression![linear regression](https://images.unsplash.com/photo-1539651044670-315229da9d2f?ixlib=rb-1.2.1&ixid=MnwxMjA3fDB8MHxwaG90by1wYWdlfHx8fGVufDB8fHx8&auto=format&fit=crop&w=1170&q=80)photo by [<NAME>](https://unsplash.com/@edanco) at [Unsplash](https://unsplash.com/) Introductionthis is the companion colab code for my tutorial in medium about Least square regression in Python with genomic data, please read the article for better understanding of the code. In the article I discuss about the theory, the hyperparameters, how to efficiently use the techniques. Here I am providing just generic information about the algorithm and the code. this is a colab notebook, you can test in colab (all the code has been tested and you do not need to set your machine)[Medium article about Least squares regression math](https://)note: the articles will be publish soonfor other stories check my medium profile:[<NAME>](https://salvatore-raieli.medium.com/)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Least regression with only one variableLet's start with one variable, then we will do with more variables for the input.<jupyter_code>n= 20
noise = 3
df = pd.DataFrame()
df["x1"] = np.linspace(-10, 10, num = n) 
df["x0"] = np.ones(n)
w = np.array([2,3])
X = df.to_numpy() # to make matrix multiplication easier
y =  w.dot(X.transpose())
y = np.random.normal(y,noise)
y
#plotting the line and the point
plt.plot(np.linspace(-10, 10, num = n), w.dot(X.transpose()))
sns.scatterplot(df["x1"], y)
# Least Square
#I am showing here all the passages
# w = (X^T * X)^-1 (X^T*t)

# X^T
X_transposte = X.transpose()

# (X^T * X)
XT_dot_X = X_transposte.dot(X)

# (X^T * X)^-1
XT_dot_X_inv = (np.linalg.inv(XT_dot_X))

# (X^T*t)
XT_dot_t = X_transposte.dot(t)

# w = (X^T * X)^-1 (X^T*t)
w = XT_dot_X_inv.dot(XT_dot_t)

# y model
y = X.dot(   w)<jupyter_output><empty_output><jupyter_text>with more than one X variablewe are generating our data, starting from X we get y, we are choosing some arbitray weight to generate the values of y. later we will try to predict using the least squares method.notice we are using some arbitrary weights, for generating the data. the Y (true values) is generating multiplying X for the weights and adding some noise (mimicking a real dataset). You can play around with the weights and to observe how similar are the weight returning by the model. Notice, the weights are similar but not perfectly the same (since we have add some noice).<jupyter_code>def dataset_gen(w = None, n = 100, noise = 0 ):
  """ generate random data for least square regression"""
  x1 = np.linspace(-10, 10, num = n) #we create random number 
  x2 = np.random.normal(x1,5) #we have to add noise otherwise we cannot invert the matrix later
  X  = np.array([np.ones(len(x1)), x1, x2]) #we are generating X0, or the  bias
  _w = np.array(w)
  y =  _w.dot(X) #we generate a y
  y = np.random.normal(y,noise) #we add some noise
  return X, y 

#we are writing a function that is doing the passage we seen before
#instead to do separatly we use a simple function
def least_square(X,y):
  """ perform least square regression"""                     
  Xt = X.dot(y)
  XX = np.linalg.inv(X.dot(X.T))
  XX = np.round(XX,4)
  Xt = np.round(Xt,2)    
  w  = XX.dot(Xt)
  w  = np.round(w,2)
   
  return w
#generate data
w = [10,2,5]
noise = 10
X, y = dataset_gen(w = w, n = 100, noise = noise )
#let's use our function and check the weight obtained
w_ls = least_square(X,y)
print(w_ls)
#generate the prediction
y_pred  = w_ls.dot(X)     
# Root of Mean Sum of Squared Error                   
RMSE = np.sqrt(((y -y_pred )**2).sum()/len(y))
RMSE
w = np.array(w)
w_ls = np.array(w_ls)
plt.plot(X[1],y,'.', label='x1 vs t')
x1 = np.linspace(-10, 10, num = 100)
X_ = np.array([np.ones(len(x1)), x1, x1])
plt.plot(x1, w.dot(X_) ,'-' , label='Generated data')
plt.plot(x1, w_ls.dot(X_) ,'--', label='Least Squares result')
plt.legend()<jupyter_output><empty_output><jupyter_text>An example with scikit-learnI will show you here a simple example with the boston dataset and using the least square method. Linear regressor estimator in Scikit-learn is using ordinary least square (so, under the hood is the same procedure, but of course you do not have to write least square function with scikit-learn). [linear regression in scikit-learn](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)the boston dataset variables:The data contains the following columns:* 'crim': per capita crime rate by town.*'zn': proportion of residential land zoned for lots over 25,000 sq.ft.*'indus': proportion of non-retail business acres per town.*'chas':Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).*'nox': nitrogen oxides concentration (parts per 10 million).*'rm': average number of rooms per dwelling.*'age': proportion of owner-occupied units built prior to 1940.*'dis': weighted mean of distances to five Boston employment centres.*'rad': index of accessibility to radial highways.*'tax': full-value property-tax rate per $10,000.*'ptratio': pupil-teacher ratio by town*'black': 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.*'lstat': lower status of the population (percent).*'medv': median value of owner-occupied homes in $$1000swe will predict the 'medv' variable.<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn import metrics
from yellowbrick.regressor import AlphaSelection, PredictionError, ResidualsPlot
%matplotlib inline
dataset = "https://raw.githubusercontent.com/SalvatoreRa/tutorial/main/datasets/Boston.csv"
df = pd.read_csv(dataset)
df = df.iloc[:,1:] #the first column is not useful (is an ID columns)
df.describe()
#prepare X and Y
X = df[['crim', 'zn', 'indus', 'chas', 'nox', 'rm', 'age', 'dis', 'rad', 'tax',
       'ptratio', 'black', 'lstat']]
y = df['medv'] #this is our target variable
#stratify
X_train, X_test, y_train, y_test = train_test_split(X , y, 
                                                        test_size = 0.2, #test size would be 20%
                                                        random_state = 42, #42 is for good luck!
) 
#noticed we have not stratified the data as in classification tasks
#stratify is for balancing the categories in the sampling, but here the variable is continuous
#scaling the white
scaler = StandardScaler()
scaler.fit(X_train)
X_train = scaler.transform(X_train)
X_test = scaler.transform(X_test) #this to avoid information leakage
model = LinearRegression() #initialize...
model.fit(X_train,y_train) #train the model...
y_pred = model.predict(X_test) #predict! 
#evaluate the model: on the evaluation metric I will go in deep in the next tutorial
print('Mean Absolute error (MAE) :', metrics.mean_absolute_error(y_test, y_pred))
print('Mean Squared Error (MSE):', metrics.mean_squared_error(y_test, y_pred))
print('Root Mean Squared Error (MSE):', np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
#visualize
plt.scatter(y_test, y_pred)
plt.xlabel('Y dataset')
plt.ylabel('LS prediction')<jupyter_output>Mean Absolute error (MAE) : 3.1890919658878483
Mean Squared Error (MSE): 24.29111947497352
Root Mean Squared Error (MSE): 4.928602182665337<jupyter_text>Bonus: nice visualization with yellowbricksThe Yellowbrick library is a python library which allows you diagnostic your model, conduct visualization and many more things. [Yellowbrick](https://www.scikit-yb.org/en/latest/)we will plot:* **Residuals Plot** which is showing the difference between real values and the predicted values. The model is also plotting the R^2 for train and test set.* **Prediction error plot** this plot is showing the true values against the prediction, the 45° line represent the optimal line, more our fitted line is far, less the model is accurate<jupyter_code># residual plot
model = LinearRegression()
visualizer = ResidualsPlot(model) #we use yellowbrick wrapper

visualizer.fit(X_train, y_train)  # we fit the model
visualizer.score(X_test, y_test)  # test set evaluation
g = visualizer.poof() #plotting!
#similarly we use the wrapper
#prediction error plot
model = LinearRegression()
visualizer = PredictionError(model)

visualizer.fit(X_train, y_train)  
visualizer.score(X_test, y_test)  
g = visualizer.poof()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10823909531502424, 'fraction_numerical': 0.022155550426955918, 'mean_word_length': 4.642578125, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 8, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '150041', 'license': 'MIT', 'n_tokens_mistral': 2990, 'n_tokens_neox': 2783, 'n_words': 1008, 'path': 'machine learning/Least_squares_regression.ipynb', 'repo_name': 'SalvatoreRa/tutorial', 'size': '8680'}