我们从Python开源项目中，提取了以下15个代码示例，用于说明如何使用scipy.special.iv()。
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']
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 _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
