How to plot scipy probability distribution from pymc3?

Questions
1

I have been looking into the bayesian interface from input data. I created the following bayesian interface model by using input data from pdf using the code from the question in the link (PDF python code),

def IQR_outlier_detection(data):
    ##### Detecting Outliers (IQR):
    ## Quantiles:
    S_list_O = np.array(data)
    PDF_Q1 = np.quantile(S_list_O, 0.25)
    PDF_Q3 = np.quantile(S_list_O, 0.75)
    print('25% of PDF =', PDF_Q1)
    print('75% of PDF =', PDF_Q3)

    IQR = PDF_Q3 - PDF_Q1
    print('IQR of PDF =', IQR)

    Lower_Quantile = PDF_Q1 - (1.5 * IQR)
    Upper_Quantile = PDF_Q3 + (1.5 * IQR)
    print('Lower limit of IQR of PDF =', Lower_Quantile)
    print('Upper limit of IQR of PDF =', Upper_Quantile)
    print(S_list_O)

    ## Removing Outliers:
    S_list_FLNO = []
    Quantile_Set = (Lower_Quantile, Upper_Quantile)
    for x in S_list_O.tolist():
        if x >= Quantile_Set[0] and x <= Quantile_Set[1]:
            S_list_FLNO.append(x)
    print('PDF with outliers = ', S_list_O)
    print('PDF with no outliers = ', S_list_FLNO)
    return S_list_FLNO

def Stats_outlier_detection(data):
    ##### Detecting Outliers (Statistics):
    ## Mean and STD of data:
    S_list_O = np.array(data)
    print(S_list_O)
    mean_S1 = np.mean(S_list_O)
    std_S1 = np.std(S_list_O)
    print('Mean of PDF =', mean_S1)
    print('STD of PDF =', std_S1)

    ## Removing Outliers:
    S_list_NO = [x for x in S_list_O if (x > mean_S1 - (2 * std_S1))]
    S_list_FLNO = [x for x in S_list_NO if (x < mean_S1 + (2 * std_S1))]
    print('PDF with outliers = ', S_list_O)
    print('PDF with no outliers = ', S_list_FLNO)
    return S_list_FLNO

def mu_and_std_beta_distribution(alpha, beta):
    # # method 1:
    # mean = (alpha * 1.0) / (alpha + beta)
    # var = (alpha * beta * 1.0) / (alpha + beta + 1.0) / (alpha + beta) ** 2.0
    # std = np.sqrt(var)

    # method 2:
    mean, var = best_dist.stats(alpha, beta, moments='mv')
    std = np.sqrt(var)
    return mean, std

def Bayesian_Interface(data):
###################################################### Outliers ########################################################
    S_FLNO = IQR_outlier_detection(data)
    # S_FLNO = Stats_outlier_detection(data)
########################################################################################################################
    # Bayesian Model:
    with pm.Model() as model_S:
        # Prior Distributions for unknown model parameters:
        alpha_S = pm.HalfNormal('alpha_S', sigma=10)
        beta_S = pm.HalfNormal('beta_S', sigma=1)

        # Observed data is from a Likelihood distributions (Likelihood (sampling distribution) of observations):
        observed_data_S = pm.Beta('observed_data_S', alpha=alpha_S, beta=beta_S, observed=S_FLNO)

        # draw 5000 posterior samples
        trace_S = pm.sample(draws=5000, tune=2000, chains=3, cores=1)

        # Obtaining Posterior Predictive Sampling:
        post_pred_S = pm.sample_posterior_predictive(trace_S, samples=3000)
        print(post_pred_S['observed_data_S'].shape)
        print('\nSummary: ')
        print(pm.stats.summary(data=trace_S))
        print(pm.stats.summary(data=post_pred_S))
    return trace_S, post_pred_S

print('#######################################################################################################')
            print('Performing Bayesian Interface for PDF:')
            trace_0, post_pred_0 = Bayesian_Interface(data=np_pdf)
            print('Accessing PDF from Bayesian Interface:')
            best_dist = getattr(st, 'beta')
            mean_0, std_0 = mu_and_std_beta_distribution(alpha=np.mean(trace_0['alpha_S']), beta=np.mean(trace_0['beta_S']))
            param_0 = [np.mean(trace_0['alpha_S']), np.mean(trace_0['beta_S']), mean_0, std_0]
            print('Parameters values = ', param_0)
            pdf_0, x_axis_pdf_0, y_axis_pdf_0 = make_pdf(dist=best_dist, params=param_0, size=len(np_pdf))
            print('#######################################################################################################')

So, I get the values from the alpha and beta of the distribution from pymc3 package in order to use these values in the same code in the link. But, I feel that I am doing something wrong. How can I create a pdf from the bayesain interface model that I am using?