You could use lmfit
(https://lmfit.github.io/lmfit-py/) and define a model function, say as
import numpy as np
from lmfit import Model
import matplotlib.pyplot as plt
x = np.array((5.107, 15.593, 178.942, 351.23, 523.172, 1039.449))
y = np.array((3.57, 4.09, 9.19, 14.3, 19.41, 32.17))
def expfunc(x, scale, decay, offset):
"model exponential decay with offset"
return offset + scale*x**decay
# create model from the above model function
model = Model(expfunc)
# create parameters with initial values,
# using names from model function
params = model.make_params(offset=0, scale=1, decay=1)
# fit data 'y' to model with params
result = model.fit(y, params, x=x)
# print and plot result
print(result.fit_report())
result.plot_fit()
plt.show()
The fit report for this fit will give fitting statistics and best-fit values, uncertainties, and correlations for the parameters:
[[Model]]
Model(expfunc)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 27
# data points = 6
# variables = 3
chi-square = 0.181
reduced chi-square = 0.060
Akaike info crit = -15.015
Bayesian info crit = -15.640
[[Variables]]
offset: 3.29036599 +/- 0.200622 (6.10%) (init= 0)
scale: 0.06290220 +/- 0.008912 (14.17%) (init= 1)
decay: 0.88280026 +/- 0.020216 (2.29%) (init= 1)
[[Correlations]] (unreported correlations are < 0.100)
C(scale, decay) = -0.997
C(offset, scale) = -0.722
C(offset, decay) = 0.686
and produce a plot like this:

FWIW, this model function is simple enough that you could use an lmfit ExpressionModel
to build the model from the string of the expression, as with:
from lmfit.models import ExpressionModel
model = ExpressionModel('offset + scale*x**decay')
with all else the same as above.