Проблемы с бинарной логистической регрессией

Передо мной стоит задача: дописать класс бинарной логистической регрессии. На скрине формулы, которые нужно использовать.

Вот код, который у меня получился, но в тестах он выдает неверный ответ. Подскажите, что здесь не так.

class LogisticRegressionGD:
    '''
    A simple logistic regression for binary classification with gradient descent
    '''

    def __init__(self):
        pass

    def __extend_X(self, X):
        """
            Данный метод должен возвращать следующую матрицу:
            X_ext = [1, X], где 1 - единичный вектор
            это необходимо для того, чтобы было удобнее производить
            вычисления, т.е., вместо того, чтобы считать X@W + b
            можно было считать X_ext@W_ext 
        """
        return np.append(X, np.ones((X.shape[0], 1)), axis=1)

    def init_weights(self, input_size, output_size):
        """
            Инициализирует параметры модели
            W - матрица размерности (input_size, output_size)
            инициализируется рандомными числами из
            нормального распределения со средним 0 и стандартным отклонением 0.01
        """
        np.random.seed(42)
        self.W = np.random.normal(scale=0.01, size=(input_size, output_size))

    def get_loss(self, p, y):
        """
            Данный метод вычисляет логистическую функцию потерь
            @param p: Вероятности принадлежности к классу 1
            @param y: Истинные метки
        """
        L = lambda p, y: -sum(y*np.log(p)+(1-y)*np.log(1-p))/y.shape[0]

        return L(p, y)


    def get_prob(self, X):
        """
            Данный метод вычисляет P(y=1|X,W)
            Возможно, будет удобнее реализовать дополнительный
            метод для вычисления сигмоиды
        """
        if X.shape[1] != self.W.shape[0]:
            X = self.__extend_X(X)

        P = 1/(1 + np.exp( -np.dot(X, self.W) ))

        return P

    def get_acc(self, p, y, threshold=0.5):
        """
            Данный метод вычисляет accuracy:
            acc = \frac{1}{len(y)}\sum_{i=1}^{len(y)}{I[y_i == (p_i >= threshold)]}
        """
        predict_class = np.zeros(p.shape[0])
        is_predict_true = np.zeros(y.shape[0])

        predict_class = np.where(p >= threshold, 1, 0)

        is_predict_true = np.where(y == predict_class, 1, 0)

        accuracy = sum(is_predict_true)/y.shape[0]

        return accuracy

    def fit(self, X, y, num_epochs=500, lr=0.001):

        X = self.__extend_X(X)
        self.init_weights(X.shape[1], y.shape[1])

        accs = []
        losses = []
        for _ in range(num_epochs):
            p = self.get_prob(X)

            # W_grad = np.dot( np.transpose(X) , self.get_prob(X) - y )/y.shape[0]
            W_grad = X.T.dot(p-y)/y.shape[0]
            self.W -= lr*W_grad

            # необходимо для стабильности вычислений под логарифмом
            p = np.clip(p, 1e-10, 1 - 1e-10)

            log_loss = self.get_loss(p, y)
            losses.append(log_loss)
            acc = self.get_acc(p, y)
            accs.append(acc)

        return accs, losses

P.S. На данном этапе еще не нужно предсказывать что-либо, нужно только оптимизировать алгоритм

P.P.S. Построил графики точности и функции потерь, вдруг это чем-то поможет.