diff --git a/machine_learning/linear_regression.py b/machine_learning/linear_regression.py
index fb5d0da5e335..99a23518ddb8 100644
--- a/machine_learning/linear_regression.py
+++ b/machine_learning/linear_regression.py
@@ -1,11 +1,11 @@
 """
 Linear regression is the most basic type of regression commonly used for
-predictive analysis. The idea is pretty simple: we have a dataset and we have
+predictive analysis. The idea is pretty simple: we have a dataset, and we have
 features associated with it. Features should be chosen very cautiously
 as they determine how much our model will be able to make future predictions.
 We try to set the weight of these features, over many iterations, so that they best
-fit our dataset. In this particular code, I had used a CSGO dataset (ADR vs
-Rating). We try to best fit a line through dataset and estimate the parameters.
+fit our dataset. In this particular code, I used a CSGO dataset (ADR vs
+Rating). We try to best fit a line through the dataset and estimate the parameters.
 """
 
 # /// script
@@ -13,17 +13,19 @@
 # dependencies = [
 #     "httpx2",
 #     "numpy",
+#     "matplotlib",
 # ]
 # ///
 
 import httpx2
+import matplotlib.pyplot as plt
 import numpy as np
 
 
 def collect_dataset():
     """Collect dataset of CSGO
     The dataset contains ADR vs Rating of a Player
-    :return : dataset obtained from the link, as matrix
+    :return : dataset obtained from the link, as a matrix
     """
     response = httpx2.get(
         "https://raw.githubusercontent.com/yashLadha/The_Math_of_Intelligence/"
@@ -41,13 +43,13 @@ def collect_dataset():
 
 
 def run_steep_gradient_descent(data_x, data_y, len_data, alpha, theta):
-    """Run steep gradient descent and updates the Feature vector accordingly_
+    """Run steep gradient descent and update the Feature vector accordingly_
     :param data_x   : contains the dataset
     :param data_y   : contains the output associated with each data-entry
     :param len_data : length of the data_
     :param alpha    : Learning rate of the model
-    :param theta    : Feature vector (weight's for our model)
-    ;param return    : Updated Feature's, using
+    :param theta    : Feature vector (weights for our model)
+    ;param return    : Updated features, using
                        curr_features - alpha_ * gradient(w.r.t. feature)
     >>> import numpy as np
     >>> data_x = np.array([[1, 2], [3, 4]])
@@ -99,19 +101,24 @@ def run_linear_regression(data_x, data_y):
 
     theta = np.zeros((1, no_features))
 
+    err = []
+
     for i in range(iterations):
         theta = run_steep_gradient_descent(data_x, data_y, len_data, alpha, theta)
         error = sum_of_square_error(data_x, data_y, theta)
         print(f"At Iteration {i + 1} - Error is {error:.5f}")
 
-    return theta
+        if i % 1000 == 0:
+            print(f"At Iteration {i + 1} - Error is {error:.5f}")
+
+    return theta, err
 
 
 def mean_absolute_error(predicted_y, original_y):
     """Return sum of square error for error calculation
     :param predicted_y   : contains the output of prediction (result vector)
     :param original_y    : contains values of expected outcome
-    :return          : mean absolute error computed from given feature's
+    :return          : mean absolute error computed from given features
 
     >>> predicted_y = [3, -0.5, 2, 7]
     >>> original_y = [2.5, 0.0, 2, 8]
@@ -122,6 +129,44 @@ def mean_absolute_error(predicted_y, original_y):
     return total / len(original_y)
 
 
+# visualization
+def plot_regression(data_x, data_y, theta):
+    """
+    Plot regression line with dataset points
+    """
+
+    x = np.array(data_x[:, 1]).flatten()
+    y = np.array(data_y).flatten()
+
+    predictions = theta[0, 0] + theta[0, 1] * x
+
+    plt.scatter(x, y)
+
+    plt.plot(x, predictions)
+
+    plt.xlabel("ADR")
+    plt.ylabel("Rating")
+
+    plt.title("Linear Regression Best Fit")
+
+    plt.show()
+
+
+def plot_loss(err):
+    """
+    Plot training loss curve
+    """
+
+    plt.plot(err)
+
+    plt.xlabel("Iterations")
+    plt.ylabel("Loss")
+
+    plt.title("Training Loss Curve")
+
+    plt.show()
+
+
 def main() -> None:
     """Driver function"""
     data = collect_dataset()
@@ -130,7 +175,11 @@ def main() -> None:
     data_x = np.c_[np.ones(len_data), data[:, :-1]].astype(float)
     data_y = data[:, -1].astype(float)
 
-    theta = run_linear_regression(data_x, data_y)
+    theta, err = run_linear_regression(data_x, data_y)
+
+    plot_regression(data_x, data_y, theta)
+    plot_loss(err)
+
     len_result = theta.shape[1]
     print("Resultant Feature vector : ")
     for i in range(len_result):