中 国卫 生 统 计 2010年 lO月 第 27卷 第 5期 最 大 熵 方 法 一功 率 谱 密 度分 析 法 在 时 间序 列 资 料 研 究 中的应 用 罗同勇 周 敦金 Ayako Sumi 余 滨 陈邦华 余松 林 【提 要】 目的 介绍一种应用较为广泛 的新 的时 间序列分析方法 “最大熵方法一功率谱密度分析法 ”。方法 对 时 间序列资料在频域分析 中采用最 大熵 方法进行功率 谱分析 ，并结合在时域分析 中对非 线性 最小二乘 法采用线性化方 法进 行分析 、预 测。结果 介绍 了该方法的理论柜架 、操 作步骤 ，并通过 英 国伦敦 1948—1967年麻 疹报 告发病 的时问序列进 行了演示 ，还对此方法与传统 的频域分析法和时域分析法进行 了比较。结论 本方法 可对 时间序列 资料 的 内在特征深入 分析 ，在 了解序列 内部结构 的基础上进行拟合预测分析 ，在 时间序列资料 分析中具有广泛的应用前景 。 【关键词 】 时 间序列 谱 分析 时域分析 最大熵方法一功率谱密度分析法 目前 在 我 国 医学 领 域 中 ，对 时 间 序 列 利 用 AR、 由决 定 分形 结 构 的 系 统 部 分 (systematic part)和 包 括 MA、ARMA 及 ARIMA 等模 型进 行 时域分 析 的文献 较 多 J。而频域分析由于其方法复杂 ，计算量大 ，在医 噪音 的波 动 部 分 (fluctuating part)组 成 ，可用 方 程 式 (1)来 表 达 。 学 领域 中 的运 用 目前在 我 国少有 文 献报 道 。本 文介 绍 (f)=systematic part+fluctuating part (1) 的最 大熵 方法 -功 率谱 密度 (MEM —PSD)分析 法 ，是一 系统 部分 被认 为是 原始 序 列 (t)的决 定 性 变量 种综合 了频 域 中利 用 最 大 熵 方 法 (maximum entropy X 、， ． (t)；波 动部 分则 是 非 决 定 性变 量 部分 如 噪音 等 ， method，MEM )进行 功 率 谱 分 析 和 时 域 分 析 中 的非 线 由原 始序列 减 去决 定性 变量 序列 后 的残差 时间序 列组 性最小 二 乘 法 (non—linear least squares method，LSM ) 成 。方程 (1)对 类似 于 图 1或 具有 分 形 结构 的时 间序 进行模 型 分 析 线 性 化 开 展 拟 合 的分 析 方 法 。本 方 列预 测分 析是 必须 的 。 法 克服 了常规 谱 分析 的弱点 和 非线 形最 小二乘 法拟 合 因此 ，时 间序 列 分 析 的关 键 就 变 成 了对 决 定 性 变 的困难 ，通过 对 时间序 列 的不 同频 率 的周 期波 动 进 行 量 的估计 。理论 上 ，时 间序 列 的 X (t)可 通 过 非 线 性 分析 ，得 出序列 的决定性 变 量特 征 ，然 后再 通过 曲线 外 最小 二乘 法 (LSM)对 决定 性 变 量 外延 得 到 。所 以 ，对 延法 得 出相应 预测 结果 。 如 图 1所 示 的决定 性 变量 (f)，可 由结合 了正 弦 和 MEM -PSD 的理论 背 景 余 弦 的方 程 得 出 ，如方 程 (2)： Np 医学 生物 学 和其 它学 科 分 支 一 样 ，均不 断 需要 有 A (t) = Olo +∑ {ansin(2~rf．t) + n = 1 新 的分析 方法来 克服 传统 分析 方 法 的缺陷 对时 间序列 b COS(2~f．O } (2) 的内在特征 即 分形结 构 进 行 深 人 分 析 ，从 而 在 深刻 了 在方程 (2)中 ，未 知参 数 包 括 、 。、OL 和 b (n= 解序列 的内部 结构 的基 础 上 ，进行 拟合 预测 分析 。 I，2，…， )， (=1／t ，t 为周 期 )为 rt次周期 波 动 的 图 I所示 时 间序列 ，为英 国伦敦 1948～1967年 麻 频 率 ，a 和 b 分别 为 第 ，z次周 期 波 动 的振 幅 ，，z 为 构 疹 报告 发病数 据 。通 过 肉眼 观察 即 能发 现 ，序 列 周 成 周期 的 总 数 ，a 为 代 表 时 间 序 列 平 均 值 的一 个 常 期 性 的呈 现 出 自相 似 性 ，我 们 将 这 种 自相 似性 称 为 分 量 。 形结 构 。 600O 4000 2000 0 1948 195l l954 l957 l960 l963 l966 图 1 英 国伦敦 1948—1967年麻疹周报告发病情况 1．对决定性变量的分离 对 于如图 l所示 时间序列 ，我们 可认为其序列值 1．武汉市疾病预防控制中 ~(430015) 2．日本札幌医科 大学 医学 院卫生系 3．华 中 科 技 大 学 同济 医学 院流 行 病 与 卫 生 统 计 教 研 室 (430030) 使用 方 程 (2)进 行 线 性最 小 二 乘 法拟 合 时必 定 是 非线 性 的 。为得 到这 些 参 数 的最 佳 优化 值 ，必 须 对 非 线性 数据 (方 程 )进 行线 性 化 。在 本 研 究 中 ，通 过 利用 最大 熵方 法估 计 的周 期 以达 到 线 性 化 ，然后 将 最 优最 ／j~--乘法 拟合 (optimum least squares fitting，LSF)得 到 的曲线延伸到时间序列 的预测范围内，就得到 (t) 的预 测值 。 2． 的确定 (最大熵方法．功率谱密度分析法 ) 值 是 通 过 最 大 熵 方 法 的 功 率 谱 密 度 (power spectral density，PSD)P( (f：频率 )的谱峰位置来确 定 的 。在抽 样 间隔 △f均 等 的条 件下 ，(功 率谱 密 度 )P △通讯作者 ：罗 同勇 ，E．mail：whnic@msn．corn (，)通过方程 (3)计算。 

· 478 · Chinese Joumal of Health Statistics．Oct2010．Vo1．27．No．5 P(，)：— — 『1+∑Tm,kexp卜 i2~rfkAt]I‘ — 一 ： 0 (3) 其 中 P 值为第 m序 次 的波 动周期 时 滤波 器 的输 出功率 ，Ym 为 当 = ，，… ( 为最 优 滤波 器 序 次 时相应 的 滤,k 0 1 m 波器 系 数 。 P 和 y ． 能通过 利 用 Burg氏 过 程解 Yule—Walker方程求得 。 3．最优 序次 值的确定 为 精确估 计 MEM—PSD，我们 需 要 一 个 合 适 的标 3．决定 。在 MEM 估 计 的周 期 的基础 上 ，指定 构成 分析期 间内基 本 变量 的周期结构 的基 础模 型。然 后，分析 l0个主要的 MEM 估计周期对所分析的最佳 最小 二乘法 配合 曲线及预 测范 围 的贡 献 。可用 以下 三 个 步骤 来完 成这一 分析 。 (1)用不 同 的 值 ，通 过 把 这 1O个 样 式 一 个 接 一 个地 加到 LSF曲线 上 以使 谱 峰频 率 的 功 率 最 大化 的方法 ，计算出分析范围内的 LSF曲线。 (2)将每一个 值计算出来 的 LSF曲线延伸至 准来确定 最优 序次 值 。在本 研究 中，除 了三个 传统 预 测 区域 。 的标 准 方 法 ：最 终 预 测 误 差 标 准 (the final prediction— (3)通 过对 LSF曲线 的评 价 ，确 定 最 优 化 最 优筛 error，VPE)、Akaike信 息理 论 标 准 (Akaike’S informa— 选 序次 M。 tion theroretical criterion，AIC)和 自回归 转换 函数标 准 (the autoregressive transfer function criterion，CAT) 外 7'引，本研究还使用 了一个称为“特征相关时间(the characteristic correlation time，CCT)”的新标 准方法 ，详 细计 算 了不 同波 动周期 数 (，，z)的 MEM—PSD峰值 ]。 4． 。、 和 b 的确定 (最 小二乘 法分 析 ) 使 用周 期 函数 即方程 (2)进 行 LSF计 算 时 ，除 非 给定这 些周期 波 动数 ，否则 将 不 可避 免 的成 为非 线 性 LSM。但 在本研 究 中 ，由于利 用 了 MEM 估 计 的周 期 数 ( )，非 线性 LSM 被线性 化 了。因此 ，除 了 Ⅳn外 ， 方程 (2)中参数 ％ 、 和 b (n=1，2，3，… ， )的最 优 值 均可 利用方 程 (2)，并采 用 MEM 求 得 的周 期数 ，结 合使用 优化 LSF计 算而得 到 。 5．确定 主要波 动周期 样式 从 残差时 间序列 (X (t)=XA(t)一x (t))的标 准 差 的趋 势 ，很 容易 地估 计 出不 同周 期样 式 对 基 本 变 量 的贡献 。从 而 ，能够 确定 主要波 动周期 NP的值 。 6．有关 “MemCalc”程序 限 制最 大 熵 方 法一功 率谱 密度 分 析方 法 广 泛应 用 的一种重 要 因素 ，就 是 以上 计 算 过程 任 务 繁重 。但 日 4．确定基 础模 型。在 贡献 率 的基础 上 ，我 们 可 以 安全 的选 择 的最 优化值 和基 础模 型 。 5．通 过 和 MEM 中估 计 的 ，利 用 方 程 2，求 出 0、 和 b 6．通 过延 伸 曲线 得到预 测值 。 时间序列 数据 预 处 理 后 的 修正后 的时间序 列数据 数 据缺 失 离 群 点 平均采样 间 对 数转换 高峰值 移 动平均 子序 列 的分 里 ●____________________________-_________●____●_●__--一 I 周期性结构的确定 l 厂————] 甄 佩两在 ! 盎 最 小二乘法 拟合 — ———— —— T—— ———一 最 佳拟合 曲线的确定 鎏裹篓塞II 溅簋 量的确定II 孝 测分析 预测 曲线I酸 差 时间序 本 学者 Tanaka Y等将 以上 频域 及 时 域 的分 析 过 程集 图2 最大熵方法．功率谱密度(MEM—PSD) 合 到名 为“MemCalc”计 算机 程 序 ，使 以上 分 析 步 骤变 得 简单 、便 捷 ¨ 。 分析法计算过程示意 图 MemCalc在 时 间序 列 中的应 用举例 操 作步 骤 我们 以英 国伦敦 1948～1967年 的 麻 疹报 告 发 病 本 方法 的操作 步骤并 简要 描述如 下 (图 2)： 资料 (图 1)为例 ，演示整 个分 析过程 。 1．建 立用 于分析 的时 间序 列数 据 。 当需 要 时 ，对 1．谱 分析 原始 数据进 行预 处理 ，如 平 均抽 样 间 隔 、填 补缺 损 值 、 将 1948～1957年 的数 据 作 为 分 析 区 间 ，将 1958 离群点 更正 、对数 转换 、数据集群 及移 动平 均数等 。 ～ 1967年作 为 预测 区间 (图 3)。通 过 “MemCalc”程 2．频 域分析 (MEM 谱分 析 )：为 了确 定预 测 滤 波 序计算 分析 区间 的 MEM—PSD。计 算 结 果 见 图 4。其 器序次 的最 优值 ，通 过最大 熵方法 ，利用 FPE、AIC、 中主要的波动周期(月)包括 ：24、12、8、46等。 CAT及 CCT等方 法 ，计 算 出每个 频率 的功 率 谱 密 的 2．确定 主要 波动周期 谱峰值 ，然后利用方程(3)计算出 值 的 MEM—PSD。 为得到 优化 的 LSF曲线 ，必须 确定 主要 波动 周 期 通过各 频率 的功率 谱密 度 ，确 定波 动 周 期 为峰 频 率 的 模 式 。我们分 别计算 了 1～10个波 动周期 模式 的 LSF 倒数 ，从 PSD 的总体 趋势得 到波动 的性质 。 曲线 ，并 计算在 每种波 动模 式下 的残差 值 ，见 图 5。从 

中 国卫 生 统 计 2010年 lO月 第 27卷 第 5期 图可见 ，当取 前 3个 波 动 周 期 后 ，标 准 差下 降较 少 ，故 (FFr)，随着研究 的不断 深入 ，以上各 方 法均 暴露 出 了 在本 次分 析 中 ，确定 此 时间序 列 的 主要 波 动 周期 为 24 局 限性 ¨…。 月 、12月和 8个月 ，然后进行 LSF曲线拟合。 20世纪 60年代 ，Burg在 分 析 地 震 信 号 时 提 出最 大熵 估计 理论 ，该理 论 克 服 了传 统谱 分 析 所 固 有 的分 妊 嚣 ． ． ．△ ．。 辨率 不高 、频 率 漏 泄 等 缺 点 ，使 谱 分 析 进 人 一 个 新 阶 段 ，使 之在 工程 、物理 学 、天文 学 、海洋 学 、气 象学 、医学 等领域 得 到广泛 运用 。 2．与常用 时域分 析 方法 的 比较 在 时 间序列 的时 域 分 析 中 ，常用 方 法包 括广 泛应 图 3 英 国伦敦 1948—1967年麻疹周报告发病情 况 用 的 ARMA、ARIMA、直 接 拟 合 技 术 等 方 法 。 与 AR— 0 ∞ l ∞ ‘ jL ． LIl^^． A II ．^ ． 一 MA 方法 比较 ，Ohtomo对 Barnes等 用 ARMA 拟 合 太 阳黑子 的分 析进行 再分 析 ，他认 为 ARMA未 能 准确 的 拟合 太 阳黑子 出现 规律 ¨。。 。与直 接拟合 技术 分析 ， 由于其 没有进 行谱 分析 ，除 非给定 周 期波 动 ，则 其拟合 过程 为非 线性 的 。 3．由于“MemCalc”程序 的发 明 ，使 MEM—PSD 分 图 4 原始数据分析区间的不同波动频 率(1／月 ) 析法 能够 得 到广 泛 利 用 。 目前 国外 已有 大 量 文献 ，利 下 的 MEM—PSD 分 布 用 MEM 方 法对传 染 病 的波 动规 律 进 行 分 析 ，从 而 对 700 600 500 200 lO0 0 图 5 不同波动周期数时残差时间序列 的标准差 ． 圆 搽 一 预测数值 l I 时 同 (年 ) 图 6 在 3个主要波动周期数时预测值与实际值 的拟合情 况 3．曲线 拟合 将确 定 的 3个 周 期数 作 为 主 要 的 波动 周 期 ，并进 未来 发病 进行预 测 。 Time Series Analysis Based by Spectral Power Distdbufion (PSD)·Maximum Entropy Method(M EM ) Luo Tongyong， Zhou Dunjin，Ayako Sumi，et a1．Wuhan Centersfo，Diseases Pre— vention and Control(430015)，Hubei 【Abstract】 Objective To introduce fl new time series analy- sis．“MEM ．PSD”．M ethods W hich is a linearized version of the nonlin— ear least squares method for fitting analysis in the time domain，combined with the maximum entropy method (MEM )for spectral analysis in the fre— quency domain．Results Theoretical framewo~ and procedures have been introduced，measles data in London have been applied tO demonstrate the whole analysis procedures，an d the comparison with the conventional methods also has been di~ussed．Conclusion Th e specialty of the meth． od is to an alyze the deterministic characteristic of time series data in fre— quency domain，it appears a bm ad application prospects． 【Key words】 Time series；Spectral analysis；Time domain 行 预测分 析 ，见 图 6。为便 于 比较 ，我们 将 实 际 观察 数 analysis；M EM -PSD 用虚 线表 示 ，将 预测 值 用 实 线表 示 ，可 见 ，预 测 值 较好 的拟合 了从 1958～1967年 实际 的发 病趋 势 。 参 考 文 献 讨 论 中国 卫 生 统 计 ，2009，26(1)：64． 1．吕效 国 ，刘 凯 峰 ，王 占君 ．用 白 回 归 模 型 的 预 选 方 法 预 测 门 诊 人 数 ． 1．与传统 谱分 析 的 比较 2．ARIMA模 型对 中国人 口死亡 率预测 的研究．中国卫生统 计 ，2008， 25(6)：630-631． 谱 分 析 即通 过计 算整个 时 间序列 频率 域 内 的功率 3．许筱红 ，金小林．自回归数学模 型对疟 疾疫情 预测的研究．中国卫生 分布 ，对 序 列 的分 形 结 构 (即 自相似 性 )进 行 分 析 ，从 统计 ，2006。23(3)：228-229． 而揭示 出被 研究 序列 中包 含 的周 期 性 ，加 深 对 序 列 的 理解 。从 2O世纪 以来 ，谱分 析方 法取 得 了几次 重大 的 4．Sumi A．，Ohtomo N．，Tanaka Y ，et a1．Prediction an alysis for measles epidemics．The Japan society of Applied Physics，2003，42．761 1-7620． 5．Tanak a Y ，Ohtomo N，Temc~ S．Nonlinear time series an alysis 2．The 突破 ，目前 常用 的谱分 析 方法 包 括 ：Schuster周期 图法 construction of a data an alysis system “M emCalc”．Bull Fac Eng Hok— (Schuster’S periodogram)、Wiener—Khintchine法则 、傅 kaido Univ，1992，158：43-55． 立叶 转 换 (Fourier transform)… 及 快 速 傅 立 叶 转 换 (下转第 484页 ) 

· 484 · 中 国卫 生 统计 2010年 1O月第 27卷 第 5期 解决 这一 问题 ，可 能 有赖 于检 测 技 术 的进 一 步 改 善 。 binomial sampfing；Delta method；Fisher information；Simulation 例 如 ，如果样 品不仅 可 以制成 混合样 本进行 检测 ，也可 以单独 进行 检测 (比如 说 ，人 或 动 物 的 血液 样 本 )，则 采用混 合样本方 法 还 可 以进 一 步 降低 成 本 。例 如 ，当 总体率 为 0．01时 ，与逐 一 检 测方 案 相 比 ，混合 样 本 方 参 考 文 献 1．姜庆五，陈启 明．流行 病 学 方 法与 模 型．上 海 ：复 旦 大学 出 版社 ， 2007，3945． 2．何书元．概率论 与数理统计．北京 ：北京 大学 出版社 ，2006，221， 法可 以减少 约 80％ 的工 作 量 ，率 越 低 ，可 以节 省 的 工 3．孙振球．医学统 计学．第 2版．北京 ：高等教育出版社 ，2006，122． 作量越 大 ¨ 。 4．顾卫东．混合样本方法 估测媒介感染 率．中国寄 生虫与 寄生虫病杂 A study on interval estim ation of rates using sequential bino- mial sampling and pooled samples Sun Qingwen，Zhang Luo— man，Yu Feifei，et a1．Department of Mathematics and P sics， Second Military Medical University(200433)，Shanghai 志，1998，16(1)：29·33． 5．Gu W ，Lampman R，Novak RJ．Assessment ofarbovirus vec tor infec tion rates using variable size pooling．Mcd Veterin Entomol，2004，1 8(2)： 200．2()4． 6．W alter SD ，Hildreth SW ，Beaty BJ．Estimation of infection rates in pop— ulations of organ isms using pools of variable size．American Journal of Epidemiology，1980，112(1)：124-128． 【Abstract】 Objective To obtain interval estimation of lower 7．Sacks JM ，Bolin SR，Crowder SV．Preyalence estimation from pooled rates using sequential binomial sampling and pooled samples for decreasing samples．Am J Vet Res，1989，50(2)：205-206． the numbers of test．M ethods 95％ approximate confident intervals are 8．俞潇潇 ，刘沛．混合检验 总体率可信 区间估计方 法．中国卫生统计 ， provided by using the central limit theorem and Delta method or Fisher in- 2oo7，24(1)：74．75． formation，an d several pooling strategies are examined from the points of 9．于振凡 ，张尧庭 ，于欣丽 ，等．数 据的统计 处理和解 释．第 2版．中国 the correctness rate，precision，and the numbers of test by theoretical an aly- 统 计 出版 社 ，2006，104． sis an d simulation using software M ATLAB7．1．R esd ~ W hile the pa- 10．李宝月 ，金欢 ，罗剑锋 ，等．负二项分 布抽样 中的患病率无偏估计．中 rameter r of nb(r yP)greaterthan 4，the correctness of95％ CI provided is 国卫生统计 ，2007，24(5)：459462． accurately round about 95％ ；For given precise of 95％ CI，the mi nimal r 11．Bernard KA ，M affei JG ，Jones SA ，et a1．W est Nile virus infection in can be approximately decided by letting m = 1，an d the given precise of birds an d mosquitoes，New York State．Emerging Infec tion Diseases， 95％ CI can be retained while r an d the numbe r of individuals of pool in— 2001，7(4)：679~85． creasing within some defmite range．Conclusion The numbers of test 12．Nasci RS，Gottfiied KL，Burkhalter KL ，et a1．Comparison of Vero cell call be evidently reduced by using po ol sam plin g m eth od an d sequential bi‘ plaque assay，TaqMan reverse transcriptase polymerase chain reaction nomi al sampling，an dforgiven precise of95％ CI，thepo oling strategies call RNA assay，an d VecTest antigen assay for detection of W est Nile virus be deduced by taking the cost of sample testing an d sam ple collecting into in field·co llec ted mosq uitoes ．Journal of the American M osquito Control account． Association，2002，18(4)：294·300． 【Key words】 Interval estimation；Pooled sample；Sequential 13．魏宗舒．概率论与数理统计教程．北京 ：高等教育出版社 ，1983，82—83． (上接第 479页 ) 9．RoveUi A ．Vulpiani A．Characteristic correlation time as estimation of 6．Childer DG．Modern Spectrum Analysis。ed．1EEE Press，New York， optimum filter length in M EM an alysis．Geophys．J R Astro Soc，1983， 1978． 72：293-306． 7．Seidou T，Ohtomo N．M EM an alysis oftime series data from combus- 10．Barnes JA ，Sargent fll HH．et a1．Sunspot cycle simulation using random tion MHD plasma．JPN J Appl Phys，1985，24：1204—1211． noise．The ancient sun ed．Peopin，R．O．，Eddy，J．A．(Pergamon Press， 8．Seidou T，Ohtomo N．M EM an alysis oftime series da ta from combus New York．1979，159-163． tion MHD plasma．II．The effect of an externally-applied magnetic field 1 1．Schuster HG．Deterministic Chaos：An troduction (Physik—Verlag， on the plasma turbulent flow．Jpn J Phys，1986，25：248-252． W einheim ，1984．