Fit and plot max likelihood Beta distribution to data

gistfile1.txt

import scipy.stats
import numpy as np
import scipy.optimize

obs_data = [
0.08982035928143713
,0.06818181818181818
,0.012987012987012988
,0.05357142857142857
,0.045454545454545456
,0.05405405405405406
,0.045112781954887216
,0.0661764705882353
,0.01948051948051948
,0.01904761904761905
,0.037037037037037035
,0.02531645569620253
,0.08391608391608392
,0.0896551724137931
,0.03252032520325203
,0.06410256410256411
,0.04705882352941177
,0.02824858757062147
,0.0392156862745098
,0.06593406593406594
,0.04285714285714286
,0.04285714285714286
,0.060109289617486336
,0.03355704697986577
,0.03664921465968586
,0.0425531914893617
,0.08433734939759036
,0.07964601769911504
,0.0650887573964497
,0.05102040816326531
,0.06896551724137931
,0.03260869565217391
,0.015625
,0.038461538461538464
,0.056818181818181816
,0.125
,0.029411764705882353
,0.044117647058823525
,0.061946902654867256
,0.015503875968992248
,0.0379746835443038
,0.0125
,0.08547008547008547
,0.06707317073170732
,0.03888888888888889
,0.02247191011235955
,0.022727272727272728
,0.0
,0.029556650246305417
,0.05027932960893855
,0.08421052631578947
,0.05521472392638037
,0.047619047619047616
,0.010000000000000002
,0.047619047619047616
,0.023255813953488372
,0.05747126436781609
,0.0693069306930693
,0.08280254777070065
,0.055900621118012424
,0.09009009009009009
,0.02702702702702703
,0.04918032786885246
,0.0684931506849315
,0.01595744680851064
,0.040983606557377046
,0.026490066225165566
,0.0380952380952381
,0.04301075268817204
,0.02803738317757009
,0.06382978723404256]


def neg_likelihood(x, *args):
    log_like = 0.0
    for i in range(len(obs_data)):
        log_like += np.log(max(scipy.stats.beta.pdf(obs_data[i], x[0],x[1]), 0.000001))
    return -1.0*log_like

#x0 = np.array([12, 20])
bnds = ((1.0, 1000), (1.0, 1000))
res = scipy.optimize.differential_evolution(neg_likelihood, bnds, maxiter=25)


import matplotlib.pyplot as plt
import seaborn as sns
import os


plt.figure()
sns.set_style('whitegrid')

sns.kdeplot(np.array(obs_data), shade=True, color="#1f77b4",
            label="Observered")

plot = sns.kdeplot(np.array(np.random.beta(3.82,74.9,10000)), shade=True, color="#fb5613",
                  label="Max likelihood Beta Distribution")

plot.set(xlim=(0))
plot.set_title("Max Likelihood Beta Fit")
fig = plot.get_figure()
img_filename = "beta_fit.png"
fig.savefig(img_filename)