Polynomial Regression Not Overfitting as Expected
Good news is that your code seems to be theoretically correct.
Bad news is that the numpy.linalg.pinv has numerical stability problem when fitting polynomials of high order. The problem is driven by decreasing determinant of the AT*A matrix.
Alternative approach would be to use numpy.linalg.lstsq which tend to provide slightly better results.
However, according to my tests with your sample data (shown below), even numpy.linalg.lstsq cannot consistently fit the points exactly as soon as number of points exceeds 10.
For example for only 15 points (and polynomial of degree 14) the determinant of AT*A decreases to ca 10E-150 making the matrix almost singular.
In the figure shown below I plot the log10 of determinant as a function of average distance between x-points across 1000 experiments: the pattern is very clear, the smaller the average distance the smaller the determinant. It appears plausible: the closer the x-points the more linearly dependent become the equations.
However, in practice, we are often interested in cases when the order of polynomials is significantly lower than the number of observations. In such cases the results are more numerically stable.

import numpy as np
import matplotlib.pyplot as plt
import math
import random
plt.rcParams['figure.figsize'] = [10, 5]
#*************************************************************
# added random seed for reproducibility
np.random.seed(123456)
#*************************************************************
##############################################################
def sample(n_samples):
    x = np.random.uniform(0, 1, size=n_samples)
    y = np.sin(2*np.pi*x) + np.random.normal(0, 0.3, n_samples)
        
    return [x,y]
##############################################################
def poly_regress(data, n_param):
    #coeff = np.zeros(n_param) # Coefficient vector initialization
    # Calculate the vandermonde matrix of input x
    vand_x = np.vander(data[0],N = n_param + 1, increasing=True)
    # Computes the coefficient vector via the normal equation
    coeff_pinv = np.linalg.pinv(vand_x.T @ vand_x) @ (vand_x.T @ data[1])
    #*************************************************************
    # added alternative more precise approach 
    coeff_lstsq = np.linalg.lstsq(vand_x, data[1])[0]
    #*************************************************************
    pred_pinv = 0
    pred_lstsq = 0
    for i in range(n_param + 1):
       pred_pinv += coeff_pinv[i]*np.power(data[0], i)
       pred_lstsq += coeff_lstsq[i]*np.power(data[0], i)
    return np.array([
        np.average(np.diff(np.sort(data[0]))), np.linalg.det(vand_x.T @ vand_x),
        np.average(abs(pred_pinv - data[1])), np.average(abs(pred_lstsq - data[1]))])            
##############################################################
# evaluate results
for num_pts in [3, 5, 15]:
    iter = 1000
    res = np.zeros(iter*4)
    res.shape = (iter,4)
    poly_order = num_pts - 1
    for ii in range(iter):
        # this function returns:
        # [0] average distance between x points: np.average(np.diff(np.sort(data[0])))
        # [1] determinant of the A^T * A matrix: np.linalg.det(vand_x.T @ vand_x)
        # [2] MAE of linalg.pinv-based approach
        # [3] MAE of linalg.lstq-based approach
        res[ii] = poly_regress(sample(num_pts), poly_order)
    plt.plot(res.T[0], np.log10(abs(res.T[1])), 'x', color='tab:grey', label= &quot;log10(abs(Det))&quot;)
    plt.xlabel(&quot;avg(diff(sort(data[0])))&quot;)
    plt.legend()
    plt.title(&quot;Num data points: &quot; + str(num_pts) + &quot;, poly degree: &quot; + str(poly_order) + 
              &quot;. Calculations repeated &quot; + str(iter) + &quot; times.&quot;)
    plt.show()
    plt.plot(res.T[0], (res.T[2]), 'rd', color='tab:red', label= &quot;MAE linalg.pinv&quot;)
    plt.plot(res.T[0], (res.T[3]), 'rP', color='tab:blue', label= &quot;MAE linalg.lstsq&quot;)
    plt.xlabel(&quot;avg(diff(sort(data[0])))&quot;)
    plt.legend()
    plt.title(&quot;Num data points: &quot; + str(num_pts) + &quot;, poly degree: &quot; + str(poly_order) + 
              &quot;. Calculations repeated &quot; + str(iter) + &quot; times.&quot;)
    plt.show()