rw.py 文件源码

def fit_gamma(ti):
    mu = ti.mean()
    sigma = ti.std()
    C = mu**2/sigma**2

    def f(x):
        return np.exp(x)/(2.*np.expm1(x)/x - 1.)

    def f1(x):
        y = 2.*np.expm1(x)/x - 1.
        z = 2./x**2*(x*np.exp(x) - np.expm1(x))
        return np.exp(x)*(1 - z/y)/y

    # newton solve
    K = 2.*C # initial value
    #print "Newton iteration:"
    for i in range(10):
        dK = -(f(K) - C)/f1(K)
        K = K + dK
    #    print i, "Residual", f(K) - C, "Value K =", K
    #print

    tau = mu*(1. - np.exp(-K))/K
    return K, tau