How to estimate the parameters of an exponential truncated power law using Python?

There is an equation of exponential truncated power law in the article below:

Gonzalez, M. C., Hidalgo, C. A., & Barabasi, A. L. (2008). Understanding individual human mobility patterns. Nature, 453(7196), 779-782.

like this:

It is an exponential truncated power law. There are three parameters to be estimated: rg0, beta and K. Now we have got several users' radius of gyration(rg), and uploaded it onto Github: radius of gyrations.txt

The following codes can be used to read data and calculate P(rg):

import numpy as np
# read radius of gyration from file
rg = []
with open('/path-to-the-data/radius of gyrations.txt', 'r') as f:
    for i in f:
        rg.append(float(i.strip('\n')))
# calculate P(rg)
rg = sorted(rg, reverse=True)
rg = np.array(rg)
prg = np.arange(len(sorted_data)) / float(len(sorted_data)-1)

or you can directly get rg and prg data as the following:

rg = np.array([ 20.7863444 ,   9.40547933,   8.70934714,   8.62690145,
     7.16978087,   7.02575052,   6.45280959,   6.44755478,
     5.16630287,   5.16092884,   5.15618737,   5.05610068,
     4.87023561,   4.66753197,   4.41807645,   4.2635671 ,
     3.54454372,   2.7087178 ,   2.39016885,   1.9483156 ,
     1.78393238,   1.75432688,   1.12789787,   1.02098332,
     0.92653501,   0.32586582,   0.1514813 ,   0.09722761,
     0.        ,   0.        ])

 prg = np.array([ 0.        ,  0.03448276,  0.06896552,  0.10344828,  0.13793103,
    0.17241379,  0.20689655,  0.24137931,  0.27586207,  0.31034483,
    0.34482759,  0.37931034,  0.4137931 ,  0.44827586,  0.48275862,
    0.51724138,  0.55172414,  0.5862069 ,  0.62068966,  0.65517241,
    0.68965517,  0.72413793,  0.75862069,  0.79310345,  0.82758621,
    0.86206897,  0.89655172,  0.93103448,  0.96551724,  1.        ])

I can plot the P(r_g) and r_g using the following python script:

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(rg, prg, 'bs', alpha = 0.3)
# roughly estimated params:
# rg0=1.8, beta=0.15, K=5
plt.plot(rg, (rg+1.8)**-.15*np.exp(-rg/5))
plt.yscale('log')
plt.xscale('log')
plt.xlabel('$r_g$', fontsize = 20)
plt.ylabel('$P(r_g)$', fontsize = 20)
plt.show()

How can I use these data of rgs to estimate the three parameters above? I hope to solve it using python.

According to @Michael 's suggestion, we can solve the problem using scipy.optimize.curve_fit

def func(rg, rg0, beta, K):
    return (rg + rg0) ** (-beta) * np.exp(-rg / K)

from scipy import optimize
popt, pcov = optimize.curve_fit(func, rg, prg, p0=[1.8, 0.15, 5])
print popt
print pcov

The results are given below:

[  1.04303608e+03   3.02058550e-03   4.85784945e+00]

[[  1.38243336e+18  -6.14278286e+11  -1.14784675e+11]
 [ -6.14278286e+11   2.72951900e+05   5.10040746e+04]
 [ -1.14784675e+11   5.10040746e+04   9.53072925e+03]]

Then we can inspect the results by plotting the fitted curve.

%matplotlib inline
import matplotlib.pyplot as plt

plt.plot(rg, prg, 'bs', alpha = 0.3)
plt.plot(rg, (rg+popt[0])**-(popt[1])*np.exp(-rg/popt[2]) )
plt.yscale('log')
plt.xscale('log')
plt.xlabel('$r_g$', fontsize = 20)
plt.ylabel('$P(r_g)$', fontsize = 20)
plt.show()