def FitGene(g, ns=20):  # for quick results subsample data
    t = time.time()
    Bsearch = list(np.linspace(0.05, 0.95, 5)) + [
        1.1
    ]  # set of candidate branching points
    GPy = (Y[g].iloc[::ns].values - Y[g].iloc[::ns].values.mean())[
        :, None
    ]  # remove mean from gene expression data
    GPt = monocle["StretchedPseudotime"].values[::ns]
    globalBranching = monocle["State"].values[::ns].astype(int)
    d = BranchedGP.FitBranchingModel.FitModel(Bsearch, GPt, GPy, globalBranching)
    print(g, "BGP inference completed in %.1f seconds." % (time.time() - t))
    # plot BGP
    fig, ax = BranchedGP.VBHelperFunctions.PlotBGPFit(
        GPy, GPt, Bsearch, d, figsize=(10, 10)
    )
    # overplot data
    f, a = PlotGene(
        monocle["State"].values,
        monocle["StretchedPseudotime"].values,
        Y[g].values - Y[g].iloc[::ns].values.mean(),
        ax=ax[0],
        s=10,
        alpha=0.5,
    )
    # Calculate Bayes factor of branching vs non-branching
    bf = BranchedGP.VBHelperFunctions.CalculateBranchingEvidence(d)["logBayesFactor"]

    fig.suptitle("%s log Bayes factor of branching %.1f" % (g, bf))
    return d, fig, ax