<jupyter_start><jupyter_text>This notebook shows an application of the gradient descent algorithm, to fit a linear regression model with 2 parameters (i.e. finding the slope and the intercept of the regression line to get the best estimation).<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as plt3D
import numpy as np
import scipy.stats as stats
from IPython.display import display

%matplotlib inline<jupyter_output><empty_output><jupyter_text>The data consists in a sample of 20 observations. For each individual, we have the age, the weight and the systolic pressure. Column "age" is removed from dataframe for simplicity<jupyter_code>df=pd.read_csv("systolic_blood_press.csv")
del df['age']
display(df.head())
display(df.describe())
m = len(df) # number of observations<jupyter_output><empty_output><jupyter_text>A scatter plot between weight and systolic pressure shows a strong correlation. This can be modelized with a regression line (y= a + bx) where a and b are the parameters to find. The instruction stats.linregress finds the best regression line so that the SSR (sum of squared residuals) is the lowest possible.<jupyter_code># Fit regression model
plt.scatter(x=df.weight, y=df.systolic_press)
res=stats.linregress(df.weight, df.systolic_press)

# Plot the regression line
abline_values = [res.slope*i+res.intercept for i in df.weight]
plt.plot(df.weight, abline_values, color="red")
plt.xlabel("Weight")
plt.ylabel("Systolic Pressure")

# Values found by scipy.stats
print("Slope from scipy.stats:", res.slope)
print("Intercept from scipy.stats:", res.intercept)<jupyter_output>Slope from scipy.stats: 1.35830370747
Intercept from scipy.stats: 30.2150330117<jupyter_text>The more the line fits, the less SSR is, and vice versa. This can be expressed as a cost function. This last a and b parameters correspond to intercept and slope respectively, and outputs the corresponding SSR. This function has a minimum SSR corresponding to the best parameter estimate<jupyter_code>def cost_function(a,b):
    error_squared=0
    # iterate through the sample and sum the squares of the distance between each point to the regression line
    for row in df.itertuples():
        index, systolic_press, weight = row
        estimated_y=a + b*weight
        error_squared += np.square(systolic_press-estimated_y)
    return error_squared/len(df)
# Visualize the cost function for parameter a (intercept)
cost_x = np.arange(res.intercept-0.5, res.intercept+0.55, 0.05)
cost_y = [cost_function(i, res.slope) for i in cost_x]

plt.plot(cost_x, cost_y)
plt.plot([res.intercept, res.intercept], [29.9, 30.30])
plt.xlabel("parameter a")
plt.ylabel("SSR/m")
# Visualize the cost function for parameter b (slope)
cost_x = np.arange(res.slope-0.5, res.slope+0.55, 0.05)
cost_y = [cost_function(res.intercept, i) for i in cost_x]

plt.plot(cost_x, cost_y)
plt.plot([res.slope, res.slope], [-10, 1750])
plt.xlabel("parameter b")
plt.ylabel("SSR/m")<jupyter_output><empty_output><jupyter_text>All the point of gradient descent is to find these minimums. Because the cost function is convex, it has a unique minimum which is local and global. Thus, one could use its derivative to find its minimum. Gradient descent starts with an initial guess and improves it at each iteration, so that it tends to the value minimizing the cost function. While approaching the minimum, the slope of the cost function tends to null, and gradients become smaller and smaller (convergence).The example below suffers from high iterations and low learning rates. The following example resolves this issue by scaling the weight in a [0,1] range<jupyter_code>def gradient_descent_iter(a,b):
    m=len(df)
    max_iteration = 5000 # very high number of iteration required to converge!
    iter = 0
    while True:
        iter += 1
        # at each gradient, it iterates through the entire sample (function sum below)
        # Note that for the derivative of the intercept a, we do not multiply by the weight, because the partial
        # derivative of a+b*weight with respect to a is 1, while in respect to b it would be weight
        deriva = sum([(a + b*df.weight[i] - df.systolic_press[i]) for i in range(m)]) / m
        derivb = sum([(a + b*df.weight[i] - df.systolic_press[i])*df.weight[i] for i in range(m)]) / m
        
        # While estimating these values separately, it converges very quickly, but together much more iterations
        # are required. The parameters influence each other. Note the different learning rates (0.3 and 0.0001)
        a = a - (0.3*deriva) 
        b = b - (0.0001*derivb)
        if (iter > max_iteration):
            return [a,b]
a,b = gradient_descent_iter(0, 0)
print("Found by gradient descent:", a,b)
print("From scipy:", res.intercept, res.slope)<jupyter_output>Found by gradient descent: 30.2150408598 1.35830361655
From scipy: 30.2150330117 1.35830370747<jupyter_text>Let's normalize the weight and see if it helps. The normalized y value "norm_df.systolic_press" must not be used, to keep predictions correct. It has been scaled here for demonstration<jupyter_code>norm_df = (df-df.min())/(df.max()-df.min())
display(norm_df)
# Fit regression model with scaled weight
plt.scatter(x=norm_df.weight, y=df.systolic_press)
res=stats.linregress(norm_df.weight, df.systolic_press)

# Plot the regression line
abline_values = [res.slope*i+res.intercept for i in norm_df.weight]
plt.plot(norm_df.weight, abline_values, color="red")
plt.xlabel("Weight")
plt.ylabel("Systolic Pressure")

# Values found by scipy.stats
print("Slope from scipy.stats:", res.slope)
print("Intercept from scipy.stats:", res.intercept)<jupyter_output>Slope from scipy.stats: 51.6155408837
Intercept from scipy.stats: 119.863077704<jupyter_text>In this example, the gradient descent is applied on a scaled weight. As a result, only 200 iterations are required, with a common learning rate of 1.5 for both a and b.<jupyter_code>def gradient_descent_normalized(a,b):
    m=len(df)
    max_iteration = 200 # much better, after scaling the weight in a [0,1] range
    alpha = 1.5
    iter = 0
    while True:
        iter += 1
        # Note that for the derivative of the intercept a, we do not multiply by the weight, because the partial
        # derivative of a+b*weight with respect to a is 1, while in respect to b it would be weight
        deriva = sum([(a + b*norm_df.weight[i] - df.systolic_press[i]) for i in range(m)]) / m
        derivb = sum([(a + b*norm_df.weight[i] - df.systolic_press[i])*norm_df.weight[i] for i in range(m)]) / m
        
        # now we case use a large learning rate alpha for both parameters
        a = a - (alpha*deriva) 
        b = b - (alpha*derivb)
        if (iter > max_iteration):
            return [a,b]
a,b = gradient_descent_normalized(0, 0)
print("Found by gradient descent:", a,b)
print("From scipy:", res.intercept, res.slope)
#%matplotlib nbagg

fig = plt.figure(num=None, figsize=(16,8))
ax = fig.add_subplot(1,1,1, projection='3d')

X = np.arange(30.215-10, 30.215, 1)
Y = np.arange(1.3583-0.1, 1.3583+0.1, 0.01)

X,Y=np.meshgrid(X,Y)
Z = [cost_function(a,b) for a,b in zip(X,Y)]
ax.plot_surface(X,Y,Z, rstride=1, cstride=1, cmap=plt.cm.RdYlGn)

plt.show()
#memo
#Slope from scipy.stats: 1.35830370747
#Intercept from scipy.stats: 30.2150330117<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10004072213926972, 'fraction_numerical': 0.032034749558843494, 'mean_word_length': 4.362445414847161, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '635694', 'license': 'MIT', 'n_tokens_mistral': 2623, 'n_tokens_neox': 2410, 'n_words': 834, 'path': 'GradientDescent_LinearRegression/GradientDescent_LR_2variables.ipynb', 'repo_name': 'sabuha/machine_learning', 'size': '7367'}