Как выбрать библиотеку Python и вид нелинейной регресии для моих данных?

У меня имеются следующие данные (mentor - преподаватель, mentee - ученик):

    mentor_cnt  mentee_cnt
0   3           3
1   3           3
2   7           7
3   18          19
4   24          27
5   37          35
6   53          56
7   63          70
8   86          89
9   102         114
10  135         149
11  154         169
12  174         202
13  232         287
14  298         386
15  343         475
16  384         552
17  446         684
18  509         883
19  469         757

Задача: сделать прогноз mentor_cnt при увеличении mentee_cnt на 500 человек.

Я подготовил линейную модель, которую не считаю надежной. Кроме того, связь между mentor_cnt и mentee_cnt явно нелинейна. Ниже приведу весь код со всеми print и графиками.

Вопрос: какую модель/модуль в Python выбрать для обучения нелинейной модели?

Есть ли смысл разбить исходные данные не по месяцам, а по неделям, чтобы выборка была больше? И как в таком случае делать прогноз месяц к месяцу?

P. S. В Statistica я бы использовал Nonlinear estimations/User specified regression, least squares. Но полагаю, это далеко не лучший выбор.

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

figure1, axis = plt.subplots(1, 3)
sns.histplot(df['mentee_cnt'], kde=True, ax=axis[0])
sns.histplot(df['mentor_cnt'], kde=True, ax=axis[1])
sns.scatterplot(data=df, y='mentee_cnt', x='mentor_cnt', ax=axis[2])
plt.show()

stat_mentee, p_mentee = stats.shapiro(df['mentee_cnt'])
stat_mentor, p_mentor = stats.shapiro(df['mentor_cnt'])
print('stat_mentee=%.3f, p_mentee=%.3f\n stat_mentor=%.3f, p_mentor=%.3f' %(stat_mentee, p_mentee, stat_mentor, p_mentor))

stat_mentee=0.826, p_mentee=0.002
stat_mentor=0.868, p_mentor=0.011

figure1, axis = plt.subplots(1, 2)
sns.boxplot(data=df, y='mentee_cnt', ax=axis[0])
sns.boxplot(data=df, y='mentor_cnt', ax=axis[1])
plt.show()

x = df['mentee_cnt'].to_numpy().reshape(-1, 1)
y = df['mentor_cnt'].to_numpy()

x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
model = LinearRegression().fit(x_train, y_train)
y_pred = model.predict(x_test)

mse = mean_squared_error(y_test, y_pred)
mae = mean_absolute_error(y_test, y_pred)
print('mse: %.3f, rmse: %.3f, mae: %.3f' % (mse, np.sqrt(mse), mae))

mse: 153.695, rmse: 12.397, mae: 10.425

print("Coefficients: \n", model.coef_)
print('Independent: \n', model.intercept_)
print("r2_score: \n", r2_score(y_test, y_pred))

Coefficients: [0.60379608]
Independent: 28.257163130866132
r2_score: 0.9948915119826979

residuals = y_test - y_pred
res = stats.probplot(residuals, plot=plt)
plt.ylabel('Residuals')
plt.xlabel('Excepted normal Value')
plt.show()

stat_residuals, p_residuals = stats.shapiro(residuals)
print('stat_residuals=%.3f, p_residuals=%.3f' %(stat_residuals, p_residuals))

stat_residuals=0.864, p_residuals=0.275

plt.rcParams['figure.figsize'] = [8, 5]
sns.histplot(data=residuals, kde=True, bins=5)
plt.show()

sns.scatterplot(y=residuals, x=y_pred)
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.show()

x_final_predict = np.append(x, x[-1] + 500).reshape(-1, 1)
y_final_predict = model.predict(x_final_predict)
round(y_final_predict[-1], 0)

787.0