Find gamma distribution parameters from empirical pdf (not from underlying data)

Question

I'm trying to estimate the parameters of a gamma distribution that fits best the following arbitrary pdf.

scipy work with MLE so ss.gamma.fit work with data such as histogram.

However, I don't have the underlying data since my distribution came from analytical calculus. (implied pdf from option chain with the Breeden-Litzenberger method, if that matters)

---

What I've tried so far:

1. Percentile method : I get alpha and beta = 1/scale but I cannot reconcile the loc parameter that ss.gamma need
2. with scipy.optimize : least_squares, minimize, curve_fit on both ss.gamma.cdf and ss.gamma.pdf to find parameters but I've nan or inf error from the solver (default parameters).
3. Percentile method to find initial guess, then minimise:

    import scipy.stats as ss
    import numpy as np
        
    x = np.linspace(1000,5500)
    y = get_implied_pdf()
    
    y_cdf = np.cumsum(y)

    perc_20 = x[np.argmin(np.abs(y_cdf - .2))]
    perc_80 = x[np.argmin(np.abs(y_cdf - .8))]
    
    alpha, beta = gamma_parameters(perc_20, 0.2, perc_80, 0.8) # from preceding blog post
    
    def f(x0):
        a, loc, scale = x0
        return np.dot(y_cdf, ss.gamma(a, loc, scale).cdf(x))
    
    res = minimize(f, [alpha,0,beta])
    plt.plot(x,y, label = 'original distribution')
    plt.plot(x,ss.gamma(*params).pdf(x), label = 'fitted distribution')
    plt.legend()

    plt.show()

The fitted parameters appears to give a flat distribution.

What I ended up so far:

I've finally ended up to generate synthetic data to fit distribution on it

    import scipy.stats as ss
    import numpy as np
        
    x = np.linspace(1000,5500)
    y = get_implied_pdf()

    y_cdf = np.cumsum(y)
    
    @np.vectorize
    def gen_data(u):
        return x[np.argmin(np.abs(y_cdf - u ))]
    
    data = gen_data(ss.uniform().rvs(10000))
    
    params = ss.gamma.fit(data)
    plt.plot(x,y, label = 'original distribution')
    plt.plot(x,ss.gamma(*params).pdf(x), label = 'fitted distribution')
    plt.legend()
    plt.show()

So, at the end, I've the parameters but it seems hackish and slow. Any ways / methods to fit ss.gamma directly on the pdf ? or any way to find better init_guess for solver ?

Thanks

edit : detail in progression to solution