相空间重构python

fromoperatorimportsubimportnumpyasnpfromsklearnimportmetricsfromsklearn.neighborsimportNearestNeighborsfromtoolzimportcurrydefglobal_false_nearest_neighbors(x,lag,min_dims=1,max_dims=10,**cutoffs):Acrossarangeofembeddingdimensions$d$,embeds$x(t)$withlag$\tau$,findsallnearestneighbors,andcomputesthepercentageofneighborsthatthatremainneighborswhenanadditionaldimensionisunfolded.See[1]formoreinformation.Parameters----------x:array-likeOriginalsignal$x(t).lag:intTimelag$\tau$inunitsofthesamplingtime$h$of$x(t)$.min_dims:int,optionalThesmallestembeddingdimension$d$totest.max_dims:int,optionalThelargestembeddingdimension$d$totest.relative_distance_cutoff:float,optionalThecutofffordeterminingneighborliness,indistanceincreaserelativetotheoriginaldistancebetweenneighboringpoints.Thedefault,15,issuggestedin[1](p.41).relative_radius_cutoff:float,optionalThecutofffordeterminingneighborliness,indistanceincreaserelativetotheradiusoftheattractor.Thedefault,2,issuggestedin[1](p.42).Returns-------dims:ndarrayThetesteddimensions$d$.gfnn:ndarrayThepercentageofnearestneighborsthatarefalseneighborsateachdimension.SeeAlso--------reconstructReferences----------[1]Arbanel,H.D.(1996).*AnalysisofObservedChaoticData*(pp.40-43).NewYork:Springer.x=_vector(x)dimensions=np.arange(min_dims,max_dims+1)false_neighbor_pcts=np.array([_gfnn(x,lag,n_dims,**cutoffs)forn_dimsindimensions])returndimensions,false_neighbor_pctsdef_gfnn(x,lag,n_dims,**cutoffs):#Globalfalsenearestneighborsataparticulardimension.#Returnspercentofallnearestneighborsthatarestillneighborswhenthenextdimensionisunfolded.#Neighborsthatcan'tbeembeddedduetolackofdataarenotcountedinthedenominator.offset=lag*n_dimsis_true_neighbor=_is_true_neighbor(x,_radius(x),offset)returnnp.mean([notis_true_neighbor(indices,distance,**cutoffs)forindices,distancein_nearest_neighbors(reconstruct(x,lag,n_dims))if(indices+offsetx.size).all()])def_radius(x):#PerArbanel(p.42):#thenominal'radius'oftheattractordefinedastheRMSvalueofthedataaboutitsmean.returnnp.sqrt(((x-x.mean())**2).mean())@currydef_is_true_neighbor(x,attractor_radius,offset,indices,distance,relative_distance_cutoff=15,relative_radius_cutoff=2):distance_increase=np.abs(sub(*x[indices+offset]))return(distance_increase/distancerelative_distance_cutoffanddistance_increase/attractor_radiusrelative_radius_cutoff)def_nearest_neighbors(y):Wrapperforsklearn.neighbors.NearestNeighbors.Yieldstheindicesoftheneighboringpoints,andthedistancebetweenthem.distances,indices=NearestNeighbors(n_neighbors=2,algorithm='kd_tree').fit(y).kneighbors(y)fordistance,indexinzip(distances,indices):yieldindex,distance[1]defreconstruct(x,lag,n_dims):Phase-spacereconstruction.Givenasignal$x(t)$,dimensionality$d$,andlag$\tau$,returnthereconstructedsignal\[\mathbf{y}(t)=[x(t),x(t+\tau),\ldots,x(t+(d-1)\tau)].\]Parameters----------x:array-likeOriginalsignal$x(t)$.lag:intTimelag$\tau$inunitsofthesamplingtime$h$of$x(t)$.n_dims:intEmbeddingdimension$d$.Returns-------ndarray$\mathbf{y}(t)$asanarraywith$d$columns.x=_vector(x)iflag*(n_dims-1)=x.shape[0]//2:raiseValueError('longestlagcannotbelongerthanhalfthelengthofx(t)')lags=lag*np.arange(n_dims)returnnp.vstack(x[lag:lag-lags[-1]orNone]forlaginlags).transpose()defami(x,y=None,n_bins=10):Calculatetheaveragemutualinformationbetween$x(t)$and$y(t)$.Parameters----------x:array-likey:array-like,optional$x(t)$and$y(t)$.Ifonly`x`ispassed,itmusthavetwocolumns;thefirstcolumndefines$x(t)$andthesecond$y(t)$.n_bins:intThenumberofbinstousewhencomputingthejointhistogram.Returns-------scalarAveragemutualinformationbetween$x(t)$and$y(t)$,innats(naturallogequivalentofbits).SeeAlso--------lagged_amiReferences----------Arbanel,H.D.(1996).*AnalysisofObservedChaoticData*(p.28).NewYork:Springer.x,y=_vector_pair(x,y)ifx.shape[0]!=y.shape[0]:raiseValueError('timeseriesmusthavethesamelength')returnmetrics.mutual_info_score(None,None,contingency=np.histogram2d(x,y,bins=n_bins)[0])deflagged_ami(x,min_lag=0,max_lag=None,lag_step=1,n_bins=10):Calculatetheaveragemutualinformationbetween$x(t)$and$x(t+\tau)$,atmultiplevaluesof$\tau$.Parameters----------x:array-like$x(t)$.min_lag:int,optionalTheshortestlagtoevaluate,inunitsofthesamplingperiod$h$of$x(t)$.max_lag:int,optionalThelongestlagtoevaluate,inunitsof$h$.lag_step:int,optionalThestepbetweenlagstoevaluate,inunitsof$h$.n_bins:intThenumberofbinstousewhencomputingthejointhistograminordertocalculatemutualinformation.See|ami|.Returns-------lags:ndarrayTheevaluatedlags$\tau_i$,inunitsof$h$.amis:ndarrayTheaveragemutualinformationbetween$x(t)$and$x(t+\tau_i)$.SeeAlso--------amiifmax_lagisNone:max_lag=x.shape[0]//2lags=np.arange(min_lag,max_lag,lag_step)amis=[ami(reconstruct(x,lag,2),n_bins=n_bins)forlaginlags]returnlags,np.array(amis)def_vector_pair(a,b):a=np.squeeze(a)ifbisNone:ifa.ndim!=2ora.shape[1]!=2:raiseValueError('withoneinput,arraymusthavebe2Dwithtwocolumns')a,b=a[:,0],a[:,1]returna,np.squeeze(b)def_vector(x):x=np.squeeze(x)ifx.ndim!=1:raiseValueError('x(t)mustbea1-dimensionalsignal')returnxTEXTimpo