clustering.py 文件源码

def test_von_mises_fits():
    N = 20000
    mu = 2.9
    kappa = 5
    mu_list = np.linspace(-np.pi,np.pi,20)
    kappa_list = range(1,50)
    mu_best_fit = []
    mu_record = []
    fig,ax = pt.subplots(nrows = 2)
    def kappa_guess_func(kappa,R_e):
        return (R_e - spec.iv(1,kappa)/spec.iv(0,kappa))**2

    def calc_best_fit(theta):
        z = np.exp(1j*theta)
        z_bar = np.average(z)
        mean_theta = np.angle(z_bar)
        R_squared = np.real(z_bar * z_bar.conj())
        R_e = np.sqrt((float(N)/(float(N)-1))*(R_squared-1./float(N)))
        tmp1 = opt.fmin(kappa_guess_func,3,args=(R_e,))
        return mean_theta, tmp1[0]

    for mu in mu_list:
        kappa_best_fit = []
        for kappa in kappa_list:
            mu_record.append(mu)
            theta = np.random.vonmises(mu,kappa,N)
            mean_theta, kappa_best = calc_best_fit(theta)
            mu_best_fit.append(mean_theta)
            kappa_best_fit.append(kappa_best)
        ax[0].plot(kappa_list, kappa_best_fit,'.')
        ax[0].plot(kappa_list, kappa_list,'-')
    ax[1].plot(mu_record,mu_best_fit,'.')
    ax[1].plot(mu_record,mu_record,'-')
    fig.canvas.draw(); fig.show()