Generalized power law using Gaussian processes

Hello,
I’ve been wrestling with a GP problem for awhile now, and I’ve hit a wall.

I’m trying to fit a generalized power law that has a GP ( f(x)) for the exponent:

y = \log(Y)\\ y = a + f(x) \log(x-c)

Exponentiating would give:

Y = A \exp(x-c)^{f(x)}

but I’ve been unable to fit f(x) properly.

Here’s the model and full notebook

with pm.Model() as model:
    ℓ = pm.Gamma("ℓ", alpha=3, beta=1)
    η = pm.HalfCauchy("η", beta=1)

    cov = η ** 2 * pm.gp.cov.Matern32(1, ℓ)
    #b = pm.Normal("b", mu=2, sigma=0.3)
    #mean = pm.gp.mean.Constant(b)
    mean = pm.gp.mean.Constant(2)
    
    gp = pm.gp.Latent(mean_func=mean, cov_func=cov)
    
    a = pm.Normal("a", mu=0, sigma=10)
    c = pm.TruncatedNormal("c", mu=0, sigma=10, shape=1, upper=x.min()-1e-6)
    f = gp.prior("f", X=x)
    σ = pm.HalfCauchy("σ", beta=1)


    mu = a + f * at.log( x - c)

    y_ = pm.Normal("y", mu=mu, sigma=σ, observed=y)

And here is the result. Ideally, f(x) would fit the slope change around x=11

and here is f(x).
image