Fitting data with a custom distribution using scipy.stats

So I noticed that there is no implementation of the Skewed generalized t distribution in scipy. It would be useful for me to fit this is distribution to some data I have. Unfortunately fit doesn't seem to be working in this case for me. To explain further I have implemented it like so

import numpy as np
import pandas as pd
import scipy.stats as st
from scipy.special import beta

class sgt(st.rv_continuous):

    def _pdf(self, x, mu, sigma, lam, p, q):

        v = q ** (-1 / p) * \
            ((3 * lam ** 2 + 1) * (
                    beta(3 / p, q - 2 / p) / beta(1 / p, q)) - 4 * lam ** 2 *
             (beta(2 / p, q - 1 / p) / beta(1 / p, q)) ** 2) ** (-1 / 2)

        m = 2 * v * sigma * lam * q ** (1 / p) * beta(2 / p, q - 1 / p) / beta(
            1 / p, q)

        fx = p / (2 * v * sigma * q ** (1 / p) * beta(1 / p, q) * (
                abs(x - mu + m) ** p / (q * (v * sigma) ** p) * (
                lam * np.sign(x - mu + m) + 1) ** p + 1) ** (
                          1 / p + q))

        return fx

    def _argcheck(self, mu, sigma, lam, p, q):

        s = sigma > 0
        l = -1 < lam < 1
        p_bool = p > 0
        q_bool = q > 0

        all_bool = s & l & p_bool & q_bool

        return all_bool

This all works fine and I can generate random variables with given parameters no problem. The _argcheck is required as a simple positive params only check is not suitable.

sgt_inst = sgt(name='sgt')
vars = sgt_inst.rvs(mu=1, sigma=3, lam = -0.1, p = 2, q = 50, size = 100)

However, when I try fit these parameters I get an error

sgt_inst.fit(vars)

RuntimeWarning: invalid value encountered in subtract
numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= fatol):

and it just returns

What I find strange is that when I implement the example custom Gaussian distribution as shown in the docs, it has no problem running the fit method.

Any ideas?