Generating random numbers from best probability distributions?
Good question. Lots of options. 

Most would recommend you use CDFs instead of histograms. So convert your observed distribution to empirical CDF (in R it’s ecdf() - I dint know python). Then, plot another line that is the theoretical CDF using your best-fit parameters over the span of a to b, where a to b gives you at least the same coverage as your original data set limits. 

Alternatively you could go for a chi-square option that’ll give you observed vs predicted and can be plotted as histogram but I don’t find that better in most cases, though it might make quantitative comparison more accessible