def normal_equation(x, y):
    """
    X: matrix of features, one sample per row (without bias unit)
    y: values (continuous) corresponding to rows (samples) in X
    """
    num_samples = y.size
    x = np.insert(x, 0, np.ones(num_samples), axis=1)

    return pinv(x)*y