白藜芦醇抗衰老作用的国际研究进展与文献分析 检索分析用词 SIRT1 and resveratrol Aging 信息分析平台 http://www.gopubmed.org/web/gopubmed/1?WEB0fepw17wls3d5I0I1I00h001000j10021000300.y Term: Aging Description: The gradual irreversible changes in structure and function of an organism that occur as a result of the passage of time. Synonyms: Biological Aging, Senescence 文献分析结果如下： Top Years Publications 2010 29 2012 23 2008 14 2009 12 2011 9 2007 7 2005 4 2006 3 2003 3 2013 2 2004 1 1 2 Top Countries Publications United States 39 China 10 Japan 8 France 8 Spain 4 South Korea 4 Australia 4 Italy 4 Taiwan 4 United Kingdom 3 Netherlands 2 Thailand 2 Canada 2 Hong Kong 1 Chile 1 Switzerland 1 Greece 1 Finland 1 Poland 1 Sweden 1 1 2 1 2 3 4 Top Cities Publications Farmington 6 Boston 4 Sydney 4 Taipei 4 Barcelona 3 Philadelphia 3 Beijing 3 New York City 3 Rochester 3 Paris 3 Tianjin 2 Bethesda 2 Fukuoka 2 Tokyo 2 Kantharalak 2 Lyon 2 Seoul 2 Cambridge 2 Montréal 2 Hong Kong 1 1 2 3 4 1 2 3 4 5 Top Journals Publications Plos One 4 Biochim Biophys Acta 4 Mech Ageing Dev 3 Aging 3 Biochem Biophys Res Commun 3 Cell 3 Arterioscler Thromb Vasc Biol 3 Nature 3 Age (dordr) 2 Free Radic Biol Med 2 J Biol Chem 2 Heart Fail Rev 2 Altern Med Rev 2 Arch Biochem Biophys 2 Curr Aging Sci 2 Int J Mol Med 2 Hum Reprod 1 Cell Metab 1 Am J Physiol Heart Circ Physiol 1 Biochem Pharmacol 1 1 2 3 4 5 1 2 3 ... 33 Top Authors Publications Morris B 3 Mukherjee S 2 Das D 2 Guarente L 2 Puigserver P 2 Sinclair D 2 Liu L 1 Liu M 1 Yin Y 1 Ye X 1 Zeng M 1 Zhao Q 1 Keefe D 1 Zhou Z 1 Liu B 1 Ghosh S 1 Yang X 1 Zheng H 1 Liu X 1 Wang Z 1 1 2 3 ... 33 1 2 3 ... 76 Top Terms Publications NAD-dependent deacetylase sirtuin-1 80 Humans 78 Stilbenes 73 Animals 73 Sirtuins 69 Longevity 67 Aging 60 Proteins 50 Mice 42 Unknown term default#fulltext 41 trihydroxystilbene synthase activity 39 Genes 37 Metabolism 37 metabolic process 37 Caloric Restriction 32 Unknown term default#review 32 NAD 31 Cell Aging 24 senescence 24 Histone Deacetylases 24 1 2 3 ... 76 publications over time world map network of top authors Top 20 of Important Phrases 1. sirtuin 1 browse it with it or it without it 85 2. in vitro browse it with it or it without it 84 3. oxidative stress browse it with it or it without it 59 4. histone deacetylase browse it with it or it without it 46 5. protein kinase browse it with it or it without it 43 6. gene expression browse it with it or it without it 38 7. caloric restriction browse it with it or it without it 30 8. transcription factor browse it with it or it without it 28 9. insulin sensitivity browse it with it or it without it 27 10. insulin resistance browse it with it or it without it 25 11. endothelial cells browse it with it or it without it 25 12. cell survival browse it with it or it without it 24 13. cell death browse it with it or it without it 22 14. stem cells browse it with it or it without it 20 15. health benefits browse it with it or it without it 19 16. cell lines browse it with it or it without it 18 17. adipose tissue browse it with it or it without it 17 18. skeletal muscle browse it with it or it without it 17 19. neuroprotective effects browse it with it or it without it 17 20. animal models browse it with it or it without it 14 1. Stilbenes/pharmacology browse it with it or it without it 267 2. Sirtuins/metabolism browse it with it or it without it 263 3. Sirtuin 1/physiology browse it with it or it without it 188 4. Sirtuin 1/metabolism browse it with it or it without it 157 5. Sirtuins/genetics browse it with it or it without it 139 6. Sirtuin 1/genetics browse it with it or it without it 85 7. Antioxidants/pharmacology browse it with it or it without it 75 8. Gene Expression Regulation/drug effects browse it with it or it without it 65 9. Stilbenes/therapeutic use browse it with it or it without it 61 10. Sirtuins/antagonists inhibitors browse it with it or it without it 49 11. Signal Transduction/drug effects browse it with it or it without it 42 12. Enzyme Inhibitors/pharmacology browse it with it or it without it 39 13. Apoptosis/drug effects browse it with it or it without it 34 14. Enzyme Activation/drug effects browse it with it or it without it 31 15. Oxidative Stress/drug effects browse it with it or it without it 30 16. Reactive Oxygen Species/metabolism browse it with it or it without it 27 17. Sirtuin 1/antagonists inhibitors browse it with it or it without it 26 18. Forkhead Transcription Factors/metabolism browse it with it or it without it 25 19. Acetylation/drug effects browse it with it or it without it 23 20. Niacinamide/pharmacology browse it with it or it without it 23 数据来源 http://pubmedpro.cn/Pubmed/Statistic 白藜芦醇(resveratrol，简称Res)是非黄酮类的 多酚化合物 ，天然的白藜芦醇在很多植物中存在，是植物为了抵御病菌入侵而产生的一种抗毒型物质。研究发现白藜芦醇是某些 草药 治疗炎症、脂类代谢紊乱和心脏疾病等的有效成分。随着对白藜芦醇的深入研究，显示白藜芦醇具有抗癌、抗氧化、 抗菌 、抗炎以及预防食物过敏等作用，可广泛地应用于医药、保健品、食品、化妆品等领域。 【新药或帮助人寿命延长至150岁】如何延长人的寿命一直是许多科学家研究的目标。据媒体报道制药巨头葛兰素史克正在测试部分新药，如果成功或许能让人的预期寿命延长至150岁。主要成分是人工合成的白藜芦醇，通过促进被称为长寿基因的SIRT1组蛋白脱乙酰酶活力，起到延缓衰老作用。 http://t.cn/zjn3LoG 英国制药巨头正在研制抗老新药 有望延寿150岁 动物测试：实验白鼠寿命延长15% 尽管针对白藜芦醇的研究已经开展了有10年之久，但相关领域的科学基础方面一直存在争议。 不过，人工合成的白藜芦醇在一些病症的临床测试中已经显现出令人满意的效果，并且这些病症包括癌症、心血管疾病、心力衰竭、2型糖尿病、老年痴呆症、帕金森症、脂肪肝、白内障、骨质疏松、关节炎等各种与年龄衰老有关的病症。 参与这一研发项目的哈佛大学遗传学教授戴维·辛克莱尔说：“从根本上来说，这些药物能够治疗某种病症，但与当前药物不同之处在于，它们还可以预防20种其他病症。实际意义上来说，它们可以延缓衰老。 ” 据悉，这些人工合成 “活化剂”已经在有限范围内使用于2型糖尿病患者和牛皮癣患者的临床治疗方面。科研人员发现，这些药物对2型糖尿病患者的新陈代谢带来了好处，并且让牛皮癣患者的皮肤减少变红的症状。 而在动物测试中，体重超标的实验白鼠在使用人工合成白藜芦醇后，其运动能力甚至达到了普通白鼠的两倍，并且寿命也长了15%。 http://news.sohu.com/20130312/n368468265.shtml 白藜芦醇抗衰老机理阐明 港大目标直指抗衰老药 http://www.biodiscover.com/news/hot/research/103322.html 纽约时报披露红酒抗衰老论文造假 http://www.biodiscover.com/news/hot/research/100452.html 科学家称发现"长生不老药" 人能活到1000岁 http://society.stnn.cc/qiwen/201010/t20101009_1428943.html
我们的结果应该投哪个期刊？ 我们近年一共发现了 3 个 比当前国内外《数理统计学》教材中普遍使用的“ Fisher z 变换”更好的初等显式 函数 。 （1） 第一个 形式 特别简单，已经被《 Transactions of Tianjin University 》录用。这是EI核心期刊。 该函数的最大误差是“Fisher z 变换”的 70.7% ，累计误差为“Fisher z 变换”的 141% 。 该 Quadratic Radical Function 很简单，对改进某重要问题的计算机算法，有很好的作用。不排除将来进入《计算机算法设计与分析》之类书籍或教材的可能性。这个是不是有点 阿Q 了 ？ （2） 第二个 形式略微复杂，已经被《 Communications in Statistics-Theory and Methods 》录用。目前是4区的SCI期刊。 该函数的最大误差是“Fisher z 变换”的 21.4% ，累计误差为“Fisher z 变换”的 8.90% 。 该 Sigmoid-like 函数首次胜过“Fisher z 变换”，不仅可以替换“Fisher z 变换”取得更好的精度，还对加深了解“Fisher z 变换”有重要的启发。 现在是 第三个 初等显式函数要投稿，不知道 投哪里 ？请您指教！ 谢谢！ 第三个函数 的最大误差约是“Fisher z 变换”的 18% ，累计误差约为“Fisher z 变换”的 4.7% 。 第三个函数的形式比第二个“应该”简单些。 除了我们的工作，Yun, Beong In 的 Approximation to the cumulative normal distribution using hyperbolic tangent based functions, Journal of the Korean Mathematical Society , 2009, 46(6): 1267-1276. 是近期的他人工作，里面有近几十年有 关工作的概述。这也是个4区的SCI期刊。 ————————— 相关背景 ————————— Sir Ronald Aylmer Fisher Photograph courtesy of Professor A W F Edwards by kind permission of Joan Fisher Box http://www.galtoninstitute.org.uk/Newsletters/GINL0306/university_of_cambridge_eugenics.htm 在《 大 英百科全书 ， Encyclopaedia Britannica 》 http://www.britannica.com/EBchecked/topic/208658/Sir-Ronald-Aylmer-Fisher Sir Ronald Aylmer Fisher, byname R.A. Fisher (born February 17, 1890, London, England—died July 29, 1962, Adelaide, Australia), British statistician and geneticist who pioneered the application of statistical procedures to the design of scientific experiments. 在《 The MacTutor History of Mathematics archive 》 http://www-history.mcs.st-andrews.ac.uk/Biographies/Fisher.html Fisher z-transform 在《苏联数学百科全书》的当前网络版 词条“Correlation (in statistics)” http://www.encyclopediaofmath.org/index.php/Correlation_(in_statistics) If one usually uses the Fisher z-transform , with replaced byz according to the formula Even at relatively small values the distribution of is a good approximation to the normal distribution with mathematical expectation and variance . On this basis one can now define approximate confidence intervals for the true correlation coefficient . For the distribution of the sample correlation ratio and for tests of the linearity hypothesis for the regression, see . References H. Cramér, "Mathematical methods of statistics" , Princeton Univ. Press (1946) B.L. van der Waerden, "Mathematische Statistik" , Springer (1957) M.G. Kendall, A. Stuart, "The advanced theory of statistics" , 2. Inference andrelationship , Griffin (1979) S.A. Aivazyan, "Statistical research on dependence" , Moscow (1968) (In Russian) 相关链接： 《胜过 Fisher z 变换！（1）》 http://bbs.sciencenet.cn/home.php?mod=spaceuid=107667do=blogid=603297 《胜过 Fisher z 变换！（2）》 http://blog.sciencenet.cn/home.php?mod=spaceuid=107667do=blogid=657534
PRO sf,time,rate,lag,sf,evenly=evenly $ ,delta=delta,logdelta=logdelta,lagsample=lagsample $ ,mincouples=mincouples ;+ ; NAME: ; sf ; ; PURPOSE: ; Compute the structure function of a lightcurve, according to the ; definitions given by Nandikotkur et al., 1997, AIP ; Conf. Proc 410, 1361P for the evenly case and Paltani et ; al., 1997 AA 327, 539P for the unevenly case. ; ; CATEGORY: ; time series analysis ; ; ; CALLING SEQUENCE: ; sf,time,rate,lag,sf ; ; INPUTS: ; time : The times at which the time series was measured ; rate : the corresponding count rates ; ; ; OPTIONAL INPUTS: ; none ; ; ; KEYWORD PARAMETERS: ; delta : width for the lag rebining in the unevenly case ; (if not set, is equal to 1e-6) ; logdelta : delta in logaritmic expression ; evenly : if set, compute the estimation of the structure ; function for evenly sample data, otherwise ; compute the estimation for unevenly sample data ; lagsample : sample of time lags, if not set, all possible ; lags in the time sample are take as lag sample. ; mincouples: minimum number of couples in a bin of the ; structure function not to be ignored, if not set, ; is equal to 1. ; ; OUTPUTS: ; sf : the "structure function"-values to each time lag ; lag : sample of time lags ; ; OPTIONAL OUTPUTS: ; none ; ; ; COMMON BLOCKS: ; none ; ; ; SIDE EFFECTS: ; none ; ; ; RESTRICTIONS: ; none ; ; ; PROCEDURE: ; none ; ; EXAMPLE: ; sf,time,rate,lag,sf,delta=2. ; ; ; MODIFICATION HISTORY: ; Version 1.0, 1998.12.21, Sara Benlloch (IAAT) ; Version 1.1, 1999.01.07, lag sample in the unevenly case modificied. ; ( benlloch@astro.uni-tuebingen.de ) ; Version 1.2, 1999.01.26, new keywords; logdelta, lagsample, ; mincoupples. ;- ;; ;; lightcurve (lc) parameters ;; ;; number of points in the time series npt = n_elements(time) ;; subtract mean from data rat = rate-mean(rate) ;; ;; Structure Function ;; IF (keyword_set(evenly)) THEN BEGIN ;; for evenly-sampled discrete time series (ti=i*bint; xi) bt = (time(npt-1)-time(0))/(npt-1) ;; Estimation of the structure function: ;; sf(lag) = 1/npt(tau)SUM_{i=0,...npt-1-lag}(x(i+lag)-x(i))^2 ;; the average been taken over all measurements separated by ;; the respective lag and npt(lag) being the number of such ;; pairs ( Nandikotkur et al., 1997, AIP Conf. Proc 410, 1361P) sf = dblarr(2,npt-1) FOR lag=1L,npt-1 DO BEGIN ; ( lag = 1,...,npt-1 ) sf = total((rat -rat )^2)/(npt-lag) sf = npt-lag END ;; Lag sample lag = bt*findgen(npt-1)+bt ENDIF ELSE BEGIN ;; for unevenly time series (ti,xi),i=1,...npt with arbritary ti ;; Estimation of the structure function: ;; sf(lag,delta)=(1/npt(lag,bin))SUM_{(i,j)/R} ^2 ;; such that npt(lag,delta) is the number of couples ;; that satisfy the relationship ;; R := ( lag-delta/2 tj-ti lag+delta/2 ) ;; ( Paltani et al., 1997 AA 327, 539P ) ;; logarithmic R := | log(|tj-ti|)-log(lag) | = delta ;; time and rate couples : tj-ti and xj-xi for all ji t_couple = time-shift(time,1) r_couple = rat-shift(rat,1) time_couples = t_couple(1:npt-1) rate_couples = r_couple(1:npt-1) FOR i=2L,npt-1 DO BEGIN t_couple = time-shift(time,i) r_couple = rat-shift(rat,i) time_couples = rate_couples = ENDFOR ;; lag sample IF keyword_set(lagsample) THEN BEGIN lag = lagsample ENDIF ELSE BEGIN lag = time_couples(sort(time_couples)) lag = lag(uniq(lag,sort(lag))) ENDELSE ;; width for the lag rebining IF (n_elements(delta) EQ 0) THEN delta=1e-6 ;; structur function estimate sf = dblarr(2,n_elements(lag)) FOR t=0L,n_elements(lag)-1 DO BEGIN IF keyword_set(logdelta) THEN BEGIN rel = abs( alog10(time_couples) - alog10(lag(t)) ) relationship = where( rel LE logdelta , count ) ENDIF ELSE BEGIN relationship = where( (time_couples LT lag(t)+delta/2.) $ AND (time_couples GT lag(t)-delta/2.), count) ENDELSE IF (count NE 0) THEN BEGIN sum = 0. sum_couples = rate_couples(relationship) FOR h=0L,count-1 DO BEGIN sum = sum + sum_couples(h)^2. sf(1,t) = (1./count)*sum ENDFOR ENDIF ELSE BEGIN sf(1,t)=0. ENDELSE sf(0,t) = count ENDFOR result = where( sf(1,*) NE 0. , count ) IF count NE 0 THEN BEGIN lag = lag(result) sf = sf(*,result) ENDIF IF keyword_set(mincouples) THEN BEGIN result = where( sf(0,*) GE mincouples , count ) IF count NE 0 THEN BEGIN lag = lag(result) sf = sf(*,result) ENDIF ENDIF ENDELSE return END
PRO acf,time,rate,lag,acf,covf=covf ;+ ; NAME: ; acf ; ; PURPOSE: ; Compute the autocorrelation function of an evenly sample lightcurve ; ; ; CATEGORY: ; time series analysis ; ; ; CALLING SEQUENCE: ; acf,time,rate,lag,acf,covf ; ; ; INPUTS: ; time : the times at which the time series was measured ; rate : the corresponding count rates ; ; ; OPTIONAL INPUTS: ; none ; ; ; KEYWORD PARAMETERS: ; none ; ; ; OUTPUTS: ; lag : sample time lags ; acf : autocorrelation values to each time lag given in ; a two-dimensional array, once without ; correction-factor and once with. ; ; OPTIONAL OUTPUTS: ; covf : autocovariance values to each time lag given in ; a one-dimensional array. ; ; COMMON BLOCKS: ; none ; ; ; SIDE EFFECTS: ; none ; ; ; RESTRICTIONS: ; evenly time series ; ; ; PROCEDURE: ; The autocorrelation function is computed according to the ; approximation of evenly sample given by ; covf(lag)=(1/(npt-lag))SUM_{j=1,..npt-lag} ; acf(lag)=covf(lag)/covf(0) ; Additionally, a corection-factor given by Sutherland et ; al. 1978, Ajp 219, 1029P, is added. ; ; EXAMPLE: ; acf,time,rate,lag,acf,covf=covf ; ; ; MODIFICATION HISTORY: ; Version 1.0, 1998.12.21, Sara Belloch IAAT, Joern Wilms IAAT. ; ( benlloch@astro.uni-tuebingen.de ) ;- ;; ;; lightcurve (lc) parameters ;; ;; dimension of lc in bins npt = n_elements(time) ;; Autocorrelation is defined for series with mean zero rat = rate-mean(rate) ;; ;; Autocovariance function estimate by ;; covf(lag)=(1/(npt-lag))SUM_{j=1,..npt-lag} ;; covf = dblarr(npt) FOR lag=0,npt-1 DO BEGIN ; ( lag = 0,...,npt-1 ) covf = total(rat *rat ) / (npt-lag) END lag = (time -time ) * findgen(npt) ;; ;; Autocorrelation function ;; acf(lag)=covf(lag)/covf(0) ;; acf = dblarr(2,npt) acf(0,*) = covf / covf(0) ;; ;; correction-factors (Sutherland et al. 1978, ApJ 219, 1029P) ;; k1 = 1. / npt k2 = 1. / (npt*npt) FOR i=0,npt-1 DO BEGIN acf(1,i) = acf(0,i) + k1 - float(i)*k2 ENDFOR END
http://www.roguewave.com/Portals/0/products/imsl-numerical-libraries/fortran-library/docs/6.0/math/default.htm?turl=evesf.htm EVESF（输出的结果为本征值的绝对值按照从大到小的顺序） Computes the largest or smallest eigenvalues and the corresponding eigenvectors of a real symmetric matrix. Required Arguments NEVEC — Number of eigenvalues to be computed. (Input) A — Real symmetric matrix of order N . (Input) SMALL — Logical variable. (Input) If . TRUE ., the smallest NEVEC eigenvalues are computed. If . FALSE ., the largest NEVEC eigenvalues are computed. EVAL — Real vector of length NEVEC containing the eigenvalues of A in decreasing order of magnitude. (Output) EVEC — Real matrix of dimension N by NEVEC . (Output) The J -th eigenvector, corresponding to EVAL ( J ), is stored in the J -th column. Each vector is normalized to have Euclidean length equal to the value one. Optional Arguments N — Order of the matrix A . (Input) Default: N = SIZE ( A ,2). LDA — Leading dimension of A exactly as specified in the dimension statement in the calling program. (Input) Default: LDA = SIZE ( A ,1). LDEVEC — Leading dimension of EVEC exactly as specified in the dimension statement in the calling program. (Input) Default: LDEVEC = SIZE ( EVEC ,1). FORTRAN 90 Interface Generic: CALL EVESF (NEVEC , A , SMALL , EVAL , EVEC ) Specific: The specific interface names are S_ EVESF and D_ EVESF . FORTRAN 77 Interface Single: CALL EVESF (N , NEVEC , A , LDA , SMALL , EVAL , EVEC , LDEVEC) Double: The double precision name is DEVESF . Description Routine EVESF computes the largest or smallest eigenvalues and the corresponding eigenvectors of a real symmetric matrix. Orthogonal similarity transformations are used to reduce the matrix to an equivalent symmetric tridiagonal matrix. Then, an implicit rational QR algorithm is used to compute the eigenvalues of this tridiagonal matrix. Inverse iteration is used to compute the eigenvectors of the tridiagonal matrix. This is followed by orthogonalization of these vectors. The eigenvectors of the original matrix are computed by back transforming those of the tridiagonal matrix. The reduction routine is based on the EISPACK routine TRED2 . See Smith et al. (1976). The rational QR algorithm is called the PWK algorithm. It is given in Parlett (1980, page 169). The inverse iteration and orthogonalization computation is discussed in Hanson et al. (1990). The back transformation routine is based on the EISPACK routine TRBAK1 . Comments 1. Workspace may be explicitly provided, if desired, by use of E5ESF/DE5ESF . The reference is: CALL E5ESF (N, NEVEC, A, LDA, SMALL, EVAL, EVEC, LDEVEC, WK, IWK) The additional arguments are as follows: WK — Work array of length 9 N . IWK — Integer array of length N . 2. Informational errors Type Code 3 1 The iteration for an eigenvalue failed to converge. The best estimate will be returned. 3 2 Inverse iteration did not converge. Eigenvector is not correct for the specified eigenvalue. 3 3 The eigenvectors have lost orthogonality. Example In this example, a DATA statement is used to set A to a matrix given by Gregory and Karney (1969, page 55). The largest two eigenvalues and their eigenvectors are computed and printed. The performance index is also computed and printed. This serves as a check on the computations. For more details, see IMSL routine EPISF . USE EVESF_INT USE EPISF_INT USE UMACH_INT USE WRRRN_INT IMPLICIT NONE ! Declare variables INTEGER LDA, LDEVEC, N PARAMETER (N=4, LDA=N, LDEVEC=N) ! INTEGER NEVEC, NOUT REAL A(LDA,N), EVAL(N), EVEC(LDEVEC,N), PI LOGICAL SMALL ! ! Set values of A ! ! A = ( 5.0 4.0 1.0 1.0) ! ( 4.0 5.0 1.0 1.0) ! ( 1.0 1.0 4.0 2.0) ! ( 1.0 1.0 2.0 4.0) ! DATA A/5.0, 4.0, 1.0, 1.0, 4.0, 5.0, 1.0, 1.0, 1.0, 1.0, 4.0, 2.0, 1.0, 1.0, 2.0, 4.0/ ! ! Find eigenvalues and vectors of A NEVEC = 2 SMALL = .FALSE. CALL EVESF (NEVEC, A, SMALL, EVAL, EVEC) ! Compute performance index PI = EPISF(NEVEC,A,EVAL,EVEC) ! Print results CALL UMACH (2, NOUT) CALL WRRRN ('EVAL', EVAL, 1, NEVEC, 1) CALL WRRRN ('EVEC', EVEC, N, NEVEC, LDEVEC) WRITE (NOUT,'(/,A,F6.3)') ' Performance index = ', PI END Output EVAL 1 2 10.00 5.00 EVEC 1 2 1 0.6325 -0.3162 2 0.6325 -0.3162 3 0.3162 0.6325 4 0.3162 0.6325 Performance index = 0.031
function =eof(X,neof) % function =eof_analysis(X,neof) % Wrapper function to perform PCA of a field X % with TWO spatial dimensions. (This code will also % work with a SINGLE spatial dimension. But it might % be easier to directly call 'principal_component_analysis'.) % This function basically transforms the data matrix into % a standard 2-d data matrix, taking into account NaNs. % Input: % X: (x,y,t) or (x,t). % neof: number of EOF/PC to return % Output: % EOFs: 3-d (x,y,e) or 2-d (x,e) matrix with spatial % patterns. % PCs: 2-d matrix (t,e) with principal components (scores) % in the columns % Var: variance of each principal component % Xrecon: 3-d (x,y,t) or 2-d (x,t) matrix with reconstructed X % WITHOUT adding back the mean % Note: Xpc is the projection of X onto the corresponding EOF. % That is, Xpc(it,ie) = nsum(nsum(X(:,:,it).*Xeofs(:,:,ie))) % Use this to project a physical field onto EOF space. % If X only has 2 dimensions, assume second dimension is time. % Insert 'y' dimension to conform with rest of function which % assumes X represents multiple realizations of a 2-D field. if ndims(X)==2 =size(X); X=reshape(X, ); end % Flatten 2-d fields into single vector =size(X); Xd=reshape(X, ); % (space,time) % remove NaNs inn=find(~isnan(Xd(:,1))); Xd=Xd(inn,:); % (space, time) n=size(Xd,1); %normalize for i=1:n % Xd(i,:)=Xd(i,:)-nanmean(Xd(i,:)); % Xd(i,:)=detrend(squeeze(Xd(i,:)),'constant'); % Xd(i,:)=detrend(squeeze(Xd(i,:))); Xd(i,:)=zscore(detrend(squeeze(Xd(i,:)),'constant')); end Xd=Xd'; Xd=double(Xd); % PCA if nargout3 =principal_component_analysis(Xd,neof); else =principal_component_analysis(Xd,neof); end % Xpcs (t,e) % Reshape EOFs to physical space Xeofs=repmat(NaN, ); Xeofs(inn,:)=EOFs; Xeofs=squeeze(reshape(Xeofs, )); Xrecon=repmat(NaN, ); for ie=1:neof for it=1:nt Xrecon(:,:,it)=Xrecon(:,:,it)+Xpcs(it,ie)*Xeofs(:,:,ie); end end if nargout3 XR=XR'; % (x,t) Xrecon=repmat(NaN, ); Xrecon(inn,:)=XR; Xrecon=squeeze(reshape(Xrecon, )); end function Xpcs=eof_project(Xanom,Xeofs,numd) % function Xpcs=eof_analysis(Xanom,Xeofs,numd) % Function to project physical field (1 or 2 spatial % dimensions) onto EOF. % Input: % Xanom: (x,y,t) or (x,t) - physical ANOMALY field % Xeofs: 3-d (x,y,e) or 2-d (x,e) matrix with spatial % patterns. % numd: if Xanom is (x,t), MUST set ndims=1 % Output: % Xpcs: 2-d matrix (t,e) with projection coefficients % (principal components) in the columns. Number % of projection coefficients returned is equal % to the number of EOFs passed in the input argument. % Samar Khatiwala (spk@ldeo.columbia.edu) if nargin3 % default is 2-d numd=2; end if ndims(Xanom)==2 numd==1 % data is (x,t) Xpcs=Xeofs'*Xanom; Xpcs=Xpcs'; else % nt=size(Xanom,3) % neofs=size(Xeofs,3) % Xpcs=repmat(NaN, ); % for it=1:nt % Xpcs(it,:)=squeeze(nansum(nansum(Xeofs.*repmat(Xanom(:,:,it), ))))'; % end % for ie=1:neofs % Xpcs(:,ie)=squeeze(nansum(nansum(Xanom.*repmat(Xeofs(:,:,ie), )))); % end % Flatten 2-d fields into single vector =size(Xanom); neofs=size(Xeofs,3); Xanom=reshape(Xanom, ); % (x,t) Xeofs=reshape(Xeofs, ); % (x,e) % remove NaNs inn=find(~isnan(Xanom(:,1))); Xanom=Xanom(inn,:); % (x,t) Xeofs=Xeofs(inn,:); % (x,e); Xpcs=Xeofs'*Xanom; Xpcs=Xpcs'; end function =principal_component_analysis(X,neof) % function =principal_component_analysis(X,neof) % Function to do a principal component analysis of % data matrix X. % Input: % X: (t,x) each row corrsponds to a sample, each column % is a variable. (Each column is a time series of a % variable.) % neof: number of EOF/PC to return % Output: % EOFs: (x,e) matrix with EOFs (loadings) in the columns % PCs: (t,e) matrix with principal components (scores) in the columns % Var: variance of each principal component % Xrecon: (t,x) reconstructed X (WITHOUT adding back the mean) % To reconstruct: Xrecon = PCs*EOFs' % Notes: (1) This routine will subtract off the mean of each % variable (column) before performing PCA. % (2) sum(var(X)) = sum(Var) = sum(diag(S)^2/(m-1)) % Samar Khatiwala (spk@ldeo.columbia.edu) if strcmp(class(X),'single') disp('WARNING: Converting input matrix X to class DOUBLE') end % Center X by subtracting off column means = size(X); %X = X - repmat(mean(X,1),m,1); r = min(m-1,n); % max possible rank of X % SVD if nargin 2 =svds(X,r); else =svds(X,min(r,neof)); end % EOFs: (x,e) % U: (t,e) % Determine the EOF coefficients PCs=U*S; % PCs=X*EOFs (t,e) % compute variance of each PC % Modified EV=diag(S).^2/sum( diag(S).^2 ); Var = EV*100; % Note: X = U*S*EOFs' % EOFs are eigenvectors of X'*X = (m-1)*cov(X) % sig^2 (=diag(S)^2) are eigenvalues of X'*X % So tr(X'*X) = sum(sig_i^2) = (m-1)*(total variance of X) if nargout3 Xrecon = PCs*EOFs'; % (t,x) end
function Y = runmean(X, m, dim, modestr) ; % RUNMEAN - Very fast running mean (aka moving average) filter % For vectors, Y = RUNMEAN(X,M) computes a running mean (also known as % moving average) on the elements of the vector X. It uses a window of % 2*M+1 datapoints. M an positive integer defining (half) the size of the % window. In pseudo code: % Y(i) = sum(X(j)) / (2*M+1), for j = (i-M):(i+M), and i=1:length(X) % % For matrices, Y = RUNMEAN(X,M) or RUNMEAN(X,M, ,1) % % - 1.33 2 3 4 4.67 % runmean( ,1,'mean') % % - 2 2 3 4 4 % runmean( ,1,1) % dimension 1 is larger than 2*(M=1)+1 ... % % - 2 4 6 8 10 % runmean(ones(10,7),3,2,'zero') ; % along columns, using mode 'zero' % runmean(repmat( ,5,1),2,2) ; % % - all NaN result % A = rand(10,10) ; A(2,7) = NaN ; % runmean(A,3,2) ; % % - column 7 is all NaN % runmean(1:2:10,100) % mean % % - 5 5 5 5 5 % % This is an incredibly fast implementation of a running mean, since % execution time does not depend on the size of the window. % % See also MEAN, FILTER % for Matlab R13 % version 3.0 (sep 2006) % Jos van der Geest % email: jos@jasen.nl % History: % 1.0 (2003) created, after a snippet from Peter Acklam (?) % 1.1 (feb 2006) made suitable for the File Exchange (extended help and % documentation) % 1.2 (feb 2006) added a warning when the window size is too big % 1.3 (feb 2006) improved help section % 2.0 (sep 2006) working across a dimension of a matrix. % 3.0 (sep 2006) several treatments of the edges. % Acknowledgements: (sep 2006) Thanks to Markus Hahn for the idea of % working in multi-dimensions and the way to treat edges. error(nargchk(2,4,nargin)) ; if ~isnumeric(m) || (numel(m) ~= 1) || (m 0) || fix(m) ~= m, error('The window size (M) should be a positive integer') ; end if nargin == 2, dim = ; else modestr = 'edge' ; end end modestr = lower(modestr) ; % check mode specifier if ~ismember(modestr,{'edge','zero','mean'}), error('Unknown mode') ; end szX = size(X) ; if isempty(dim), dim = min(find(szX1)) ; end if m == 0 || dim ndims(X), % easy Y = X ; else mm = 2*m+1 ; if mm = szX(dim), % if the window is larger than X, average all sz2 = ones(size(szX)) ; sz2(dim) = szX(dim) ; Y = repmat(mean(X,dim),sz2) ; else % here starts the real stuff % shift dimensions so that the desired dimensions comes first = shiftdim(X, dim-1); szX = size(X) ; % make the rest of the dimensions columns, so we have a 2D matrix % (suggested of Markus Hahn) X = reshape(X,szX(1), ; % the cumsum trick (by Peter Acklam ?) Y = cumsum(Y,1) ; Y = (Y(mm+1:end,:)-Y(1:end-mm,:)) ./ mm ; % reshape into original size Y = reshape(Y,szX) ; % and re-shift the dimensions Y = shiftdim(Y,ndims(Y)-nshifts) ; end end % ===================== % CODE OF VERSION 1.3 % ===================== % function Y = runmean(X,m) ; % % RUNMEAN - Very fast running mean filter for vectors % % Y = RUNMEAN(X,M) computes a running mean on vector X using a window of % % 2*M+1 datapoints. X is a vector, and M an positive integer defining % % (half) the size of the window. In pseudo code: % % Y(i) = sum(X(j)) / (2*M+1), for j = (i-M):(i+M), and i=1:length(X) % % % % If the total window size (2M+1) is larger than the length of the vector, the overall % % average is returned. % % % % Example: % % runmean(1:10,1) % - % % % % % % This is an incredibly fast implementation of a running average, since % % execution time does not depend on the size of the window. % % % % X should not contains NaNs (a NaN will result in a all NaN result) % % At both ends the values of Y can be inaccurate, as the first and last % % values of X are used multiple times. % % % % See also MEAN % % % for Matlab R13 % % version 1.3 (feb 2006) % % Jos van der Geest % % email: jos@jasen.nl % % % History: % % 1.0 (2003) created, after a snippet from Peter Acklam (?) % % 1.1 (feb 2006) made suitable for the File Exchange (extended help and % % documentation) % % 1.2 (feb 2006) added a warning when the window size is too big % % 1.3 (feb 2006) improved help section % % error(nargchk(2,2,nargin)) ; % % sz = size(X) ; % % if numel(sz) ~= 2 || (min(sz) ~= 1), % error('X should be a vector') ; % end % % if any(isnan(X)), % error('NaNs cannot be dealt with') ; % end % % if ~isnumeric(m) || (numel(m) ~= 1) || (m 0) || fix(m) ~= m, % error('The window size (M) should be a positive integer') ; % elseif m == 0, % Y = X ; % return ; % end % % mm = 2*m+1 ; % % if mm = prod(sz), % % if the window is larger than X, average all % warning('Window size is larger than the length of the vector.') % Y = repmat(mean(X),sz) ; % else % % the cumsum trick ... % Y = ; % Y = ; % Y = (Y(mm+1:end)-Y(1:end-mm)) / mm ; % Y = reshape(Y,sz) ; % end
function p=gkdeb(x,p) % GKDEB Gaussian Kernel Density Estimation with Bounded Support % % Usage: % p = gkdeb(d) returns an estmate of pdf of the given random data d in p, % where p.pdf and p.cdf are the pdf and cdf vectors estimated at % p.x locations, respectively and p.h is the bandwidth used for % the estimation. % p = gkdeb(d,p) specifies optional parameters for the estimation: % p.h - bandwidth % p.x - locations to make estimation % p.uB - upper bound % p.lB - lower bound. % p.alpha - to calculate inverse cdfs at p.alpha locations % % Without output, gkdeb(d) and gkdeb(d,p) will disply the pdf and cdf % (cumulative distribution function) plot. % % See also: hist, histc, ksdensity, ecdf, cdfplot, ecdfhist % Example 1: Normal distribution %{ gkdeb(randn(1e4,1)); %} % Example 2: Uniform distribution %{ clear p p.uB=1; p.lB=0; gkdeb(rand(1e3,1),p); %} % Example 3: Exponential distribution %{ clear p p.lB=0; gkdeb(-log(1-rand(1,1000)),p); %} % Example 4: Rayleigh distribution %{ clear p p.lB=0; gkdeb(sqrt(randn(1,1000).^2 + randn(1,1000).^2),p); %} % V3.2 by Yi Cao at Cranfield University on 7th April 2010 % % Check input and output error(nargchk(1,2,nargin)); error(nargoutchk(0,1,nargout)); n=length(x); % Default parameters if nargin2 N=100; h=median(abs(x-median(x)))/0.6745*(4/3/n)^0.2; xmax=max(x); xmin=min(x); xmax=xmax+3*h; xmin=xmin-3*h; dx=(xmax-xmin)/(N-1); p.x=xmin+(0:N-1)*dx; p.pdf=zeros(1,N); p.cdf=zeros(1,N); p.h=h; dxdz=ones(size(p.x)); z=p.x; else =checkp(x,p); N=numel(p.x); h=p.h; end % Gaussian kernel function kerf=@(z)exp(-z.*z/2); ckerf=@(z)(1+erf(z/sqrt(2)))/2; nh=n*h*sqrt(2*pi); for k=1:N p.pdf(k)=sum(kerf((p.x(k)-x)/h)); p.cdf(k)=sum(ckerf((p.x(k)-x)/h)); end p.x=z; p.pdf=p.pdf.*dxdz/nh; dx= ; p.cdf=p.cdf/n; % p.cdf=cumsum(p.pdf.*dx); if isfield(p,'alpha') n=numel(p.alpha); p.icdf=p.alpha; for k=1:n alpha=p.alpha(k); ix=find(p.cdfalpha,1)-1; x1=p.x(ix); x2=p.x(ix+1); F1=p.cdf(ix); F2=p.cdf(ix+1); p.icdf(k)=x1+(alpha-F1)*(x2-x1)/(F2-F1); end end % Plot if ~nargout subplot(211) plot(p.x,p.pdf,'linewidth',2) grid % set(gca,'ylim', ) ylabel('f(x)') title('Estimated Probability Density Function'); subplot(212) plot(p.x,p.cdf,'linewidth',2) ylabel('F(x)') title('Cumulative Distribution Function') xlabel('x') grid meanx = sum(p.x.*p.pdf.*dx); varx = sum((p.x-meanx).^2.*p.pdf.*dx); text(min(p.x),0.6,sprintf('mean(x) = %g\n var(x) = %g\n',meanx,varx)); if isfield(p,'alpha') numel(p.alpha)==1 text(min(p.x),0.85,sprintf('icdf at %g = %g',p.alpha,p.icdf)); end end function =checkp(x,p) n=numel(x); %check structure p if ~isstruct(p) error('p is not a structure.'); end if ~isfield(p,'uB') p.uB=Inf; end if ~isfield(p,'lB') p.lB=-Inf; end if p.lB-Inf || p.uBInf =bounded(x,p); else if ~isfield(p,'h') p.h=median(abs(x-median(x)))/0.6745*(4/3/n)^0.2; end error(varchk(eps, inf, p.h, 'Bandwidth, p.h is not positive.')); if ~isfield(p,'x') N=100; xmax=max(x); xmin=min(x); xmax=xmax+3*p.h; xmin=xmin-3*p.h; dx=(xmax-xmin)/(N-1); p.x=xmin+(0:N-1)*dx; end dxdz=ones(N,1); z=p.x; end p.pdf=zeros(size(p.x)); p.cdf=zeros(size(p.x)); function =bounded(x,p) if p.lB==-Inf dx=@(t)1./(p.uB-t); y=@(t)-log(p.uB-t); zf=@(t)(p.uB-exp(-t)); elseif p.uB==Inf dx=@(t)1./(t-p.lB); y=@(t)log(t-p.lB); zf=@(t)exp(t)+p.lB; else dx=@(t)(p.uB-p.lB)./(t-p.lB)./(p.uB-t); y=@(t)log((t-p.lB)./(p.uB-t)); zf=@(t)(exp(t)*p.uB+p.lB)./(exp(t)+1); end x=y(x); n=numel(x); if ~isfield(p,'h') p.h=median(abs(x-median(x)))/0.6745*(4/3/n)^0.2; end h=p.h; if ~isfield(p,'x') N=100; xmax=max(x); xmin=min(x); xmax=xmax+3*h; xmin=xmin-3*h; p.x=xmin+(0:N-1)*(xmax-xmin)/(N-1); z=zf(p.x); else z=p.x; p.x=y(p.x); end dxdz=dx(z); function msg=varchk(low,high,n,msg) % check if variable n is not between low and high, returns msg, otherwise % empty matrix if n=low n=high msg=[]; end
标签: function
