def nena(locs, info, callback=None):
bin_centers, dnfl_ = next_frame_neighbor_distance_histogram(locs, callback)
def func(d, a, s, ac, dc, sc):
f = a * (d / s**2) * _np.exp(-0.5 * d**2 / s**2)
fc = ac * (d / sc**2) * _np.exp(-0.5 * (d**2 + dc**2) / sc**2) * _iv(0, d * dc / sc)
return f + fc
pdf_model = _lmfit.Model(func)
params = _lmfit.Parameters()
area = _np.trapz(dnfl_, bin_centers)
median_lp = _np.mean([_np.median(locs.lpx), _np.median(locs.lpy)])
params.add('a', value=area/2, min=0)
params.add('s', value=median_lp, min=0)
params.add('ac', value=area/2, min=0)
params.add('dc', value=2*median_lp, min=0)
params.add('sc', value=median_lp, min=0)
result = pdf_model.fit(dnfl_, params, d=bin_centers)
return result, result.best_values['s']
python类iv()的实例源码
def _log_vMF_matrix(points,means,K,n_jobs=1):
"""
This method computes the log of the density of probability of a von Mises Fischer law. Each line
corresponds to a point from points.
:param points: an array of points (n_points,dim)
:param means: an array of k points which are the means of the clusters (n_components,dim)
:param cov: an array of k arrays which are the covariance matrices (n_components,dim,dim)
:return: an array containing the log of density of probability of a von Mises Fischer law (n_points,n_components)
"""
n_points,dim = points.shape
n_components,_ = means.shape
dim = float(dim)
log_prob = K * np.dot(points,means.T)
# Regularisation to avoid infinte terms
bessel_term = iv(dim*0.5-1,K)
idx = np.where(bessel_term==np.inf)[0]
bessel_term[idx] = np.finfo(K.dtype).max
log_C = -.5 * dim * np.log(2*np.pi) - np.log(bessel_term) + (dim/2-1) * np.log(K)
return log_C + log_prob
def _log_vMF_matrix(points,means,K,n_jobs=1):
"""
This method computes the log of the density of probability of a von Mises Fischer law. Each line
corresponds to a point from points.
:param points: an array of points (n_points,dim)
:param means: an array of k points which are the means of the clusters (n_components,dim)
:param cov: an array of k arrays which are the covariance matrices (n_components,dim,dim)
:return: an array containing the log of density of probability of a von Mises Fischer law (n_points,n_components)
"""
n_points,dim = points.shape
n_components,_ = means.shape
dim = float(dim)
log_prob = K * np.dot(points,means.T)
# Regularisation to avoid infinte terms
bessel_term = iv(dim*0.5-1,K)
idx = np.where(bessel_term==np.inf)[0]
bessel_term[idx] = np.finfo(K.dtype).max
log_C = -.5 * dim * np.log(2*np.pi) - np.log(bessel_term) + (dim/2-1) * np.log(K)
return log_C + log_prob
def _EM_GMM_expectation_step(self,):
self.probs = self.zij*0#np.ones((self.instance_array.shape[0],self.n_clusters),dtype=float)
#instance_array_c and instance_array_s are used to speed up cos(instance_array - mu) using
#the trig identity cos(a-b) = cos(a)cos(b) + sin(a)sin(b)
#this removes the need to constantly recalculate cos(a) and sin(a)
for mu_tmp, std_tmp, p_hat, cluster_ident in zip(self.mu_list,self.std_list,self.pi_hat,range(self.n_clusters)):
#norm_fac_exp = self.n_clusters*np.log(1./(2.*np.pi)) - np.sum(np.log(spec.iv(0,kappa_tmp)))
#norm_fac_exp = -self.n_dimensions*np.log(2.*np.pi) - np.sum(np.log(spec.iv(0,std_tmp)))
norm_fac_exp = self.n_dimensions*np.log(1./np.sqrt(2.*np.pi)) + np.sum(np.log(1./std_tmp))
#pt1 = kappa_tmp * (self.instance_array_c*np.cos(mu_tmp) + self.instance_array_s*np.sin(mu_tmp))
pt1 = -(self.instance_array - mu_tmp)**2/(2*(std_tmp**2))
self.probs[:,cluster_ident] = p_hat * np.exp(np.sum(pt1,axis=1) + norm_fac_exp)
prob_sum = (np.sum(self.probs,axis=1))[:,np.newaxis]
self.zij = self.probs/(prob_sum)
#Calculate the log-likelihood - note this is quite an expensive computation and not really necessary
#unless comparing different techniques and/or checking for convergence
#This is to prevent problems with log of a very small number....
#L = np.sum(zij[probs>1.e-20]*np.log(probs[probs>1.e-20]))
#L = np.sum(zij*np.log(np.clip(probs,1.e-10,1)))
#############################################################################
#####################Plotting functions#####################################
def pdf_bessel(self, tt):
a = self.K
t = tt/self.tau
return np.exp(-t)*np.sqrt(a/t)*iv(1., 2.*np.sqrt(a*t))/self.tau/np.expm1(a)
def vonMisesPDF(self, alpha, mu=0.0):
""" Probability density function of Von Mises pdf for input points
@param List[float] alpha List-like or single value of input values
@return List[float] List of probabilities for input points
"""
num = np.exp(self.kappa * np.cos(alpha - mu))
den = 2 * np.pi * iv(0, self.kappa)
return num / den
def kappa_guess_func(kappa,R_e):
return (R_e - spec.iv(1,kappa)/spec.iv(0,kappa))**2
def _EM_VMM_expectation_step_hard(mu_list, kappa_list, instance_array):
start_time = time.time()
n_clusters = len(mu_list);
n_datapoints, n_dimensions = instance_array.shape
probs = np.ones((instance_array.shape[0],n_clusters),dtype=float)
for mu_tmp, kappa_tmp, cluster_ident in zip(mu_list,kappa_list,range(n_clusters)):
#We are checking the probability of belonging to cluster_ident
probs_1 = np.product(np.exp(kappa_tmp*np.cos(instance_array-mu_tmp))/(2.*np.pi*spec.iv(0,kappa_tmp)),axis=1)
probs[:,cluster_ident] = probs_1
assignments = np.argmax(probs,axis=1)
#return assignments, L
return assignments, 0
def generate_bessel_lookup_table(self):
self.kappa_lookup = np.linspace(0,100,10000)
self.bessel_lookup_table = [spec.iv(1,kappa_lookup)/spec.iv(0,kappa_lookup), kappa_lookup]
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()
def convert_kappa_std(kappa,deg=True):
'''This converts kappa from the von Mises distribution into a
standard deviation that can be used to generate a similar normal
distribution
SH: 14June2013
'''
R_bar = spec.iv(1,kappa)/spec.iv(0,kappa)
if deg==True:
return np.sqrt(-2*np.log(R_bar))*180./np.pi
else:
return np.sqrt(-2*np.log(R_bar))
def generate_bessel_lookup_table(self):
self.kappa_lookup = np.linspace(0,100,10000)
self.bessel_lookup_table = [spec.iv(1,kappa_lookup)/spec.iv(0,kappa_lookup), kappa_lookup]
def _EM_VMM_GMM_expectation_step(self,):
#Can probably remove this and modify zij directly
self.probs = self.zij*0
#instance_array_c and instance_array_s are used to speed up cos(instance_array - mu) using
#the trig identity cos(a-b) = cos(a)cos(b) + sin(a)sin(b)
#this removes the need to constantly recalculate cos(a) and sin(a)
for mu_tmp, kappa_tmp, mean_tmp, std_tmp, p_hat, cluster_ident in zip(self.mu_list,self.kappa_list,self.mean_list, self.std_list, self.pi_hat,range(self.n_clusters)):
#norm_fac_exp = self.n_clusters*np.log(1./(2.*np.pi)) - np.sum(np.log(spec.iv(0,kappa_tmp)))
norm_fac_exp_VMM = -self.n_dimensions*np.log(2.*np.pi) - np.sum(np.log(spec.iv(0,kappa_tmp)))
norm_fac_exp_GMM = -self.n_dimensions_amps*np.log(2.*np.pi) - np.sum(np.log(std_tmp))
pt1_GMM = -(self.instance_array_amps - mean_tmp)**2/(2*(std_tmp**2))
#Note the use of instance_array_c and s is from a trig identity
#This speeds this calc up significantly because mu_tmp is small compared to instance_array_c and s
#which can be precalculated
pt1_VMM = kappa_tmp * (self.instance_array_c*np.cos(mu_tmp) + self.instance_array_s*np.sin(mu_tmp))
self.probs[:,cluster_ident] = p_hat * np.exp(np.sum(pt1_VMM,axis=1) + norm_fac_exp_VMM +
np.sum(pt1_GMM,axis=1) + norm_fac_exp_GMM)
#This is to make sure the row sum is 1
prob_sum = (np.sum(self.probs,axis=1))[:,np.newaxis]
self.zij = self.probs/(prob_sum)
#Calculate the log-likelihood - note this is quite an expensive computation and not really necessary
#unless comparing different techniques and/or checking for convergence
#This is to prevent problems with log of a very small number....
#L = np.sum(zij[probs>1.e-20]*np.log(probs[probs>1.e-20]))
#L = np.sum(zij*np.log(np.clip(probs,1.e-10,1)))
def _EM_VMM_expectation_step_soft(mu_list, kappa_list, instance_array, pi_hat, c_arr = None, s_arr = None):
n_clusters = len(mu_list);
n_datapoints, n_dimensions = instance_array.shape
probs = np.ones((instance_array.shape[0],n_clusters),dtype=float)
#c_arr and s_arr are used to speed up cos(instance_array - mu) using
#the trig identity cos(a-b) = cos(a)cos(b) + sin(a)sin(b)
#this removes the need to constantly recalculate cos(a) and sin(a)
if c_arr==None or s_arr==None:
c_arr = np.cos(instance_array)
s_arr = np.sin(instance_array)
# c_arr2 = c_arr[:,np.newaxis,:]
# s_arr2 = s_arr[:,np.newaxis,:]
# pt1 = (np.cos(mu_list))[np.newaxis,:,:]
# pt2 = (np.sin(mu_list))[np.newaxis,:,:]
# pt2 = (c_arr2*pt1 + s_arr2*pt2)
# pt2 = (np.array(kappa_list))[np.newaxis,:,:] * pt2
#print pt2.shape
for mu_tmp, kappa_tmp, p_hat, cluster_ident in zip(mu_list,kappa_list,pi_hat,range(n_clusters)):
norm_fac_exp = len(mu_list[0])*np.log(1./(2.*np.pi)) - np.sum(np.log(spec.iv(0,kappa_tmp)))
pt1 = kappa_tmp * (c_arr*np.cos(mu_tmp) + s_arr*np.sin(mu_tmp))
probs[:,cluster_ident] = p_hat * np.exp(np.sum(pt1,axis=1) + norm_fac_exp)
#old way without trig identity speed up
#probs[:,cluster_ident] = p_hat * np.exp(np.sum(kappa_tmp*np.cos(instance_array - mu_tmp),axis=1)+norm_fac_exp)
#older way including everything in exponent, and not taking log of hte constant
#probs[:,cluster_ident] = p_hat * np.product( np.exp(kappa_tmp*np.cos(instance_array-mu_tmp))/(2.*np.pi*spec.iv(0,kappa_tmp)),axis=1)
prob_sum = (np.sum(probs,axis=1))[:,np.newaxis]
zij = probs/((prob_sum))
#This was from before using np.newaxis
#zij = (probs.T/prob_sum).T
#Calculate the log-likelihood - note this is quite an expensive computation and not really necessary
#unless comparing different techniques and/or checking for convergence
valid_items = probs>1.e-20
#L = np.sum(zij[probs>1.e-20]*np.log(probs[probs>1.e-20]))
L = np.sum(zij[valid_items]*np.log(probs[valid_items]))
#L = np.sum(zij*np.log(np.clip(probs,1.e-10,1)))
return zij, L
def EM_VMM_clustering(instance_array, n_clusters = 9, n_iterations = 20, n_cpus=1, start='random'):
#This is for the new method...
instance_array_complex = np.exp(1j*instance_array)
kappa_lookup = np.linspace(0,100,10000)
bessel_lookup_table = [spec.iv(1,kappa_lookup)/spec.iv(0,kappa_lookup), kappa_lookup]
print '...'
n_dimensions = instance_array.shape[1]
iteration = 1
#First assignment step
mu_list = np.ones((n_clusters,n_dimensions),dtype=float)
kappa_list = np.ones((n_clusters,n_dimensions),dtype=float)
LL_list = []
if start=='k_means':
print 'Initialising clusters using a fast k_means run'
cluster_assignments, cluster_details = k_means_clustering(instance_array, n_clusters=n_clusters, sin_cos = 1, number_of_starts = 1)
mu_list, kappa_list, convergence_mu, convergence_kappa = _EM_VMM_maximisation_step_hard(mu_list, kappa_list, instance_array, cluster_assignments, instance_array_complex = instance_array_complex, bessel_lookup_table= bessel_lookup_table, n_cpus=n_cpus)
elif start=='EM_GMM':
print 'Initialising clusters using a EM_GMM run'
cluster_assignments, cluster_details = EM_GMM_clustering(instance_array, n_clusters=n_clusters, sin_cos = 0, number_of_starts = 1)
mu_list, kappa_list, convergence_mu, convergence_kappa = _EM_VMM_maximisation_step_hard(mu_list, kappa_list, instance_array, cluster_assignments, instance_array_complex = instance_array_complex, bessel_lookup_table = bessel_lookup_table, n_cpus=n_cpus)
else:
print 'Initialising clusters using random start points'
mu_list = np.random.rand(n_clusters,n_dimensions)*2.*np.pi - np.pi
kappa_list = np.random.rand(n_clusters,n_dimensions)*20
cluster_assignments, L = _EM_VMM_expectation_step_hard(mu_list, kappa_list,instance_array)
while np.min([np.sum(cluster_assignments==i) for i in range(len(mu_list))])<20:#(instance_array.shape[0]/n_clusters/4):
print 'recalculating initial points'
mu_list = np.random.rand(n_clusters,n_dimensions)*2.*np.pi - np.pi
kappa_list = np.random.rand(n_clusters,n_dimensions)*20
cluster_assignments = _EM_VMM_expectation_step_hard(mu_list, kappa_list,instance_array)
print cluster_assignments
convergence_record = []
converged = 0;
while (iteration<=n_iterations) and converged!=1:
start_time = time.time()
cluster_assignments, L = _EM_VMM_expectation_step_hard(mu_list, kappa_list,instance_array)
LL_list.append(L)
mu_list, kappa_list, convergence_mu, convergence_kappa = _EM_VMM_maximisation_step_hard(mu_list, kappa_list, instance_array, cluster_assignments, instance_array_complex = instance_array_complex, bessel_lookup_table = bessel_lookup_table, n_cpus=n_cpus)
print 'Time for iteration %d :%.2f, mu_convergence:%.3f, kappa_convergence:%.3f, LL: %.8e'%(iteration,time.time() - start_time,convergence_mu, convergence_kappa,L)
convergence_record.append([iteration, convergence_mu, convergence_kappa])
if convergence_mu<0.01 and convergence_kappa<0.01:
converged = 1
print 'Convergence criteria met!!'
iteration+=1
print 'AIC : %.2f'%(2*(mu_list.shape[0]*mu_list.shape[1])-2.*LL_list[-1])
cluster_details = {'EM_VMM_means':mu_list, 'EM_VMM_kappas':kappa_list, 'EM_VMM_LL':LL_list}
return cluster_assignments, cluster_details
