Кусочно-постоянная функция

from random import random
 
# Генерация тестовых данных
 
def make_data(a,b,c,sa,sb,n,xmin,xmax):
    r=[]
    for i in range(n):
        # генерируем случайную точку xmax > x > xmin
        x=random()*(xmax-xmin)+xmin
        # генерируем соотв. y
        if x<c:
            y=random()*sa*0.5+a
        else:
            y=random()*sb*0.5+b
        r.append((x,y))
    return r
 
# среднее по кластеру
 
def calc_avg(clast):
    s=sum(map(lambda pair: pair[1],clast))
    avg=s/len(clast)
    return avg
 
# вот данные, соответствующие ступеньке при c=3 a=-1 b=1
# 1000 точек, равномерно разбросанных от -5 до 5 (0.3 - разброс по y)
 
dd=sorted(make_data(-1,1, 3 ,0.3,0.3,1000,-5,5))
 
avg_prev,da_prev=None,None
 
for i in range(10,1000,10):
    clast_left=dd[0:i]
    avg_curr=calc_avg(clast_left)
    if avg_prev is None:
        avg_prev=avg_curr
    else:
        da_curr=abs(avg_curr-avg_prev)
        if da_prev is None:
            da_prev=da_curr
        else:
            if da_curr/da_prev>30: # резкое изменение 
                print(i)
            da_prev=da_curr
            avg_prev=avg_curr

import matplotlib.pyplot as plt
import numpy as np
n=5
# исходные данные
y=np.random.uniform(0, 3, size=n)
x=np.random.uniform(-10, 10, size=n)
# сортировака данных по x
idx=np.argsort(x)
x=x[idx]
y=y[idx]
 
a=[]
b=[]
c=[]
r=[]
for i in range(1, n):
    c.append(x[i])
    a.append(np.mean(y[:i]))
    b.append(np.mean(y[i:]))
    r.append(np.std(y[:i])+np.std(y[i:]))
idx=np.argmin(r)
print(a[idx], b[idx], c[idx])
plt.plot(x,y, 'ro')
stepx=np.array([x[0],c[idx],c[idx], x[-1]])
stepy=np.array([a[idx], a[idx],b[idx],b[idx]])
plt.plot(stepx, stepy, 'g-')

import numpy as np
from sklearn.metrics import mean_squared_error as mse
%matplotlib inline
import matplotlib.pyplot as plt
def solve(A):
    optimal_params = {k: 0 for k in 'abc'}
    best_mse = None
    domain = A[A[:, 0].argsort()]
 
    for i in range(1, len(domain) - 1):
        y_left, y_right = domain[:i, 1], domain[i:, 1]
 
        mean_left, mean_right = np.mean(y_left), np.mean(y_right)
 
        loss = np.mean(np.square(np.concatenate([y_left - mean_left, y_right - mean_right])))
        if (best_mse is None or loss < best_mse):
            best_mse = loss
            optimal_params['a'], optimal_params['b'], optimal_params['c'] = mean_left, mean_right, domain[i, 0]
            
    return optimal_params['a'], optimal_params['b'], optimal_params['c']
 
# для проверки решения
def create_f(a, b, c):
    return lambda x: a if x < c else b
 
def calc_loss(domain, f):
    return mse(domain[:, 1], np.array([f(x) for x in domain[:, 0]]))
 
n=100
# исходные данные
y=np.random.uniform(0, 3, size=n)
x=np.random.uniform(-10, 10, size=n)
# сортировака данных по x
idx=np.argsort(x)
x=x[idx]
y=y[idx]
 
A = np.stack([x, y]).T
 
a=[]
b=[]
c=[]
r=[]
for i in range(1, n):
    c.append(x[i])
    a.append(np.mean(y[:i]))
    b.append(np.mean(y[i:]))
    r.append(np.std(y[:i]) + np.std(y[i:])) # в сумме std для подмножеств y[:i] и y[i:] сладагемые участвуют с разными "коэфициентами": 1/i и 1/(n - i) соответственно, в отличие от mse, где все слагаемые участвуют с 1/n. Минимум суммы std не совпадает с минимумом mse (строго убедиться в этом можно, взяв производные) и, следовательно, не решает задачу.
idx=np.argmin(r)
print(a[idx], b[idx], c[idx])
a, b, c = a[idx], b[idx], c[idx]
calc_loss(A, create_f(a, b, c)) # mse для данных a, b, c
 
a, b, c = solve(A)
calc_loss(A, create_f(a, b, c)) # mse для данных a, b, c
 
plt.plot(x, y, 'ro')
stepx=np.array([x[0], c, c, x[-1]])
stepy=np.array([a, a, b, b])
plt.plot(stepx, stepy, 'g-');