I guess you need to look for a non-parametric approach, such as generalised additive models (GAM), e.g. regression splines. These method(s) are extremely flexible in fitting non-linear data. With polynomials, you still work in a parametric world. With regression splines, you are much more flexible and there is no need to worry about the parameterisation of the model.
I have not worked with regression splines in Python, but there are libraries. You can get a good overview of the method(s) in Chapter 7 of "Introduction to Statistical Learning". There also is Python code for the Labs in the ISL-book. So you can directly adapt these methods.
Here is a little R example:
# Lead packages
library(ISLR)
library(gam)
library(Metrics)
# OLS
reg1 = lm(mpg~qsec,data=mtcars)
# Generalised additive models (regression splines, 5 DF)
reg2 = gam(mpg~s(qsec,5),data=mtcars)
par(mfrow=c(1,1))
plot(reg2)
mae(mtcars$mpg, predict(reg1, newdata=mtcars))
mae(mtcars$mpg, predict(reg2, newdata=mtcars))
The simple GAM has an MAE of 3.4 while the linear (OLS) model has an MAE of 4.2. So quite an improvement.
This would be the GAM plot of the simple model above, including CI bands. Rather flexible fit with no effort.
