34、数据处理与分析中的关键方法及注意事项

数据处理与分析中的关键方法及注意事项

1. 功率谱密度与AR1过程匹配

在实际操作中，我们会调整(σ,β)的选择，使得AR1过程的功率谱密度（p.s.d.）大致与数据的功率谱密度相匹配。数据功率谱密度的95%置信区间可近似为AR1过程的置信区间，即cC95/jv( f )j2 。这里，C95约为6，它是自由度为2的卡方分布随机变量的95%置信水平。近似式s2 + 2σs2 也非常准确。

2. 自助法置信区间

在数据分析中，有许多特定用途的统计检验方法，但数据场景极为多样，很多情况并没有成熟的检验方法。当模型参数与数据的关系复杂，常规的误差传播方法难以应用时，确定置信区间就变得困难。

假设通过复杂的数据分析过程确定了一个模型参数m。如果有大量重复数据集（即不同时间进行相同测量得到的数据集），就可以通过经验方法来确定参数m的置信区间。对每个数据集进行相同分析，构建参数m的直方图，当重复数据集足够多时，直方图将近似于概率密度函数p(m)，进而可从p(m)推导出置信区间。虽然真实的重复数据集很少，但如果能从单个可用数据集中构建近似的重复数据集，该方法同样适用。

一个数据集d可看作是概率密度函数p(d)的N个实现。我们可以通过无限次复制这N个实现并混合，构建另一个概率密度函数p0(d)。当N趋于无穷大时，p0(d)趋近于p(d)，只要N足够大，p0(d)就是p(d)的合适近似。

这表明可以通过对原始数据集进行有放回的随机重采样来创建近似的重复数据集。若原始数据为列向量dorig，新数据集d可通过从dorig中随机选择元素N次来构建。由于允许重复选择，新数据集与原始数据集不同，且可以构建多个不同的重采样数据集。对每个重采样数据集进行相同分析，组装估计模型参数的直方图，这一过程称为自助法。

以确定直线拟合数据的斜率b的置信区间为例，我们已知如何确定该线性问题的置信区间，这是验证自助法结果的好方法。斜率的概率密度函数p(b)是正态分布，方差为σb2 ，95%置信区间在均值的2σb范围内。

自助法包含两个步骤，在对原始数据重采样Nr次的循环中进行：
- 第一步：将数据dorig和对应的时间torig重采样为d和t。
- 第二步：使用标准方法（这里是最小二乘法）从d和t中估计感兴趣的参数（这里是斜率b）。

以下是MATLAB和Python的实现代码：

MATLAB代码 ：

% edama_12_14: bootstrap statistics of slope
...
for p = [1:Nr] % loop over realizations
    % resample
    rowindex=unidrnd(N,N,1);
    t = torig(rowindex);
    d = dorig(rowindex);
    % straight line fit
    M=2;
    G=zeros(N,M);
    G(:,1)=1;
    G(:,2)=t;
    mest=(G'*G)\(G'*d);
    slope(p)=mest(2);
end

Python代码 ：

# edapy_12_14: bootstrap statistics of slope
...
for p in range(Nr):
    # resample
    # N random integers, each between 1 and N-1
    rowindex=np.random.randint(N,size=N);
    # random resampling of time with duplication
    t = eda_cvec(torig[rowindex,0]);
    # random resampling of data with duplication
    d = eda_cvec(dorig[rowindex,0]);
    # straight line fit
    M=2;
    # number of model parameters
    G=np.zeros((N,M));
    # data akernel setup
    G[:,0]=np.ones((N,));
    G[:,1]=t.ravel();
    GTG = np.matmul(G.T,G);
    # G-transpose times G
    mest=la.solve(GTG,np.matmul(G.T,d)); # least sq soln
    slope[p,0]=mest[1,0];
    # save slope in vector

重采样索引数组rowindex由MATLAB的unidrnd()函数和Python的np.random.randint()方法创建，返回1到N之间的随机整数。最终得到长度为Nr的斜率列向量，可据此形成直方图并转换为概率密度估计p(b)。

MATLAB代码 ：

% edama_12_14: bootstrap statistics of slope
...
Nbins = 100;
[shist, bins]=hist(slope, Nbins);
Db = bins(2)-bins(1);
pbootstrap = shist / (Db*sum(shist));

Python代码 ：

Nhist=1000;
# number of bins
slopemin = 0.2;
# staring bin
slopemax = 0.8;
# ending bin
ct, ed = np.histogram(slope,Nhist,(slopemin,slopemax)); # histogram
Nc = len(ct);
# length of counts ct
Ne = len(ed);
# length of edges ed
slopehist = eda_cvec(ct); # counts as a vector
# bin centers as a vector
centers = eda_cvec(0.5*(ed[1:Ne]+ed[0:Ne-1]));
# approximate pdf via normalized histogram
Ds = centers[1,0] - centers[0,0];
# bin sampling
norm = Ds*np.sum(slopehist);
# normalization factor
pbootstrap =slopehist/norm;
# p.d.f. from histogram

最后一个命令将直方图转换为面积为1的归一化概率密度函数pbootstrap。由于概率密度函数近似正态分布，95%置信区间可表示为b ± 2σb 。但在其他情况下，概率密度函数可能非正态，此时需要明确计算置信区间。

MATLAB代码 ：

% edama_12_14: bootstrap statistics of slope
...
% 5% and 95% confidence of bootstrap distribution
Pbootstrap = Db*cumsum(pbootstrap);
ilo = find(Pbootstrap>=0.025,1);
ihi = find(Pbootstrap>=0.975,1);
blo = bins(ilo);
bhi = bins(ihi);

Python代码 ：

edapy_12_14: bootstrap statistics of slope
...
# 5% and 95% confidence of
# cumulative probability function
Pbootstrap = eda_cvec(Ds*np.cumsum(pbootstrap));
# find 0.025 cumulative probability
ilovec = np.where(Pbootstrap>=0.025);
# np.where() returns a bunch of stuff,
# but what we want, ilo, is the first thing
ilox = ilovec[0];
ilo = ilox[0];
blo = centers[ilo];
# find 0.975 cumulative probability
ihivec = np.where(Pbootstrap>=0.975);
ihix = ihivec[0];
ihi = ihix[0]
bhi = centers[ihi];

这里使用MATLAB的cumsum()函数和Python的np.cumsum()方法将概率密度函数积分成概率分布Pbootstrap，使用MATLAB的find()函数和Python的np.where()方法找到包含95%总概率的斜率值blo和bhi，95%置信区间为blo < b < bhi。

在更实际的例子中，对于大西洋岩石数据集的方差最大因子分析，如果因子2的CaO与Na2O比率r对分析结论很重要，由于数据与r的关系复杂，使用标准方法推导置信区间不切实际，而自助法估计p(r)及比率r的置信区间则很直接。

3. 广义最小二乘法步骤

广义最小二乘法步骤如下表所示：
|步骤|操作|
|----|----|
|1|用文字描述问题，即数据与模型的关系|
|2|将问题整理成标准形式，确定数据d（长度N）和模型参数m（长度M），定义数据核G使dobs = Gm|
|3|检查数据，绘制数据图|
|4|确定数据的准确性，根据测量技术的准确性设定先验方差σd2 |
|5|用文字陈述先验信息，如模型参数接近已知值hpri、模型参数均值接近已知值、模型参数随空间和时间平滑变化等|
|6|将先验信息整理成标准形式：hpri = Hm|
|7|确定先验信息的准确性，根据先验信息的准确性设定先验方差σh2 |
|8|估计模型参数mest及其协方差Cm，mest = (FTF)-1FTfobs，Cm = (FTF)-1 ，其中F = [σd-1G; σh-1H]，fobs = [σd-1dobs; σh-1hpri]|
|9|陈述估计值及其95%置信区间，mtruei = mesti ± 2σmi ，σmi = √Cm[ii]|
|10|检查个体误差，计算dpre = Gmest，e = dobs - dpre，hpre = Hmest，ep = hpri - hpre，绘制ei与i、epi与i的图，绘制di_pre与di_obs、hi_pre与hi_pri的散点图，检查是否有异常大的误差|
|11|检查总误差ET，ET = E + Ep，E = σd-2eTe，Ep = σh-2eTpep，使用卡方检验评估零假设（ET仅因随机变化而与预期不同）的可能性|
|12|如果有两个不同模型，使用F检验评估零假设（两个模型的E值仅因随机变化而不同）的可能性|

4. 计算功率谱密度步骤

计算功率谱密度的步骤如下：
1. 计算时间参数 ：
- 时间序列长度d必须为偶数N。
- 时间采样间隔△t。
- 时间向量t = [0, △t, 2△t, ⋯, (N - 1)△t]T 。
- 时间序列持续时间T = N△t。
2. 计算频率参数 ：
- 奈奎斯特（最大）频率f max = 1/(2△t)。
- 频率采样间隔△f = f max/(N/2)。
- 非负频率数量Nf = (N/2 + 1)。
- 非负频率向量f = [0, △f, 2△f, ⋯, f max]T 。
3. 预处理时间序列d ：
- 减去均值，将dn替换为dn - d。
- 时间序列功率P = var(d) = N-1∑n=1Ndn2 。
- 应用汉明窗函数，将dn替换为wndn，wn = 0.54 - 0.46cos(2πtn/(T - △t))。
- 窗函数功率ff = N-1∑n=1Nwn2 。
4. 计算傅里叶变换ed ：ed = △t FFT(d)。
5. 保留非负频率 ：将ed截断为长度Nf。
6. 计算功率谱密度s2 ：s2n = (2/T)|ed[n]|2 。
7. 检查功率谱密度并查找峰值 ：绘制sn2 与fn的图。
8. 陈述零假设 ：频谱峰值是由零均值、方差为var(d)的不相关时间序列（或参数为β的AR1时间序列）的随机波动引起的。
9. 计算95%置信水平 ：
- 缩放常数c = ffvar(d)/(2Nf△f)。
- 不相关95%置信水平C95 = 6c。
- AR1 95%置信水平C95(fn) = 6c/|vn|2 ，|vn|2 = 1 + β2 - 2βcos(πfn/fny)。
10. 绘制置信水平并评估峰值的显著性 。

5. 变量持久性问题

在MATLAB和Python中，变量会永久存在，不仅可以在创建它们的单元格中访问，还可以在后续执行的其他单元格中访问，无论该单元格在脚本中的位置如何。这种特性有好处，例如可以将长脚本拆分成多个逻辑段，逐段调试，但也会导致一些常见错误：
- 变量在单元格中未正确定义，未报告错误，而是使用了之前执行的单元格中同名变量的值，该值在当前上下文中可能不正确。
- 定义变量的代码行意外从脚本中删除，但脚本仍能运行，因为变量在脚本的早期版本中已定义。但在关闭并重新打开MATLAB LiveScript或Jupyter Notebook后，脚本会意外失败。
- 预定义常量（如pi）的值被之前的脚本意外重置为意外值。MATLAB和Python都允许pi = 2这样的语句，将π的值重新定义为2。函数和方法也可以被重新定义，这种意外重新定义通常会生成难以诊断的错误消息。例如，在MATLAB中设置zeros = 3，然后执行x = zeros(3, 1)会生成“Index exceeds matrix dimensions”错误；在Python中设置np = 3（重新定义整个Numpy模块），然后执行x = np.zeros((3, 1))会生成“AttributeError: ‘int’ object has no attribute ‘zeros’”错误。

为避免或检测这些问题，可以在运行脚本前删除工作区中的所有变量。在MATLAB中，可在单元格开头放置clear all命令；Jupyter Notebook有类似命令%reset -f，但它会删除导入的模块和用户定义的变量，因此可将该命令放在00单元格开头，执行时会删除所有内容并重新导入所有模块。此外，还可以使用MATLAB的clear x和Python的del x命令删除单个变量，这在处理占用大量内存的大矩阵时很有用。

6. 时间格式转换

数据集通常包含以日历（年、月、日）和时钟（小时、分钟、秒）格式表示的时间，这种格式不适合数据分析，需要转换为用单个均匀递增变量t表示的格式。时间变量t的单位选择和起始时间t = 0的定义取决于数据分析的需求。

以黑岩森林数据集为例，该数据集包含12.6年的每小时样本，以数据集第一年1月1日为起始的以天为单位的时间变量是合理的选择，特别是考虑到昼夜周期很强。但在检查年度周期性时，以数据集第一年1月1日为起始的以年为单位的时间可能更合适，此时起始时间应便于识别特定测量的季节。

由于月份长度不同和闰年等特殊情况，将日历/时钟时间转换为单个变量t很复杂。MATLAB提供了时间算术函数，Python也有类似方法，可加快转换过程。它们接受日历日期（年、月、日）和时间（小时、分钟、秒），返回可用于计算时间差的信息。例如，计算2008年2月11日03:04:00和2008年2月11日03:04:01之间的秒数：

MATLAB代码 ：

% edapy_13_02: time arithmetic
...
DT = 86400*(datenum(2008,2,11,4,4,1)- ...
datenum(2008,2,11,4,4,0));
fprintf('DT: %.2f\n',DT);

综上所述，在数据处理与分析中，我们需要掌握功率谱密度匹配、自助法置信区间估计、广义最小二乘法、变量管理和时间格式转换等关键方法和技术，同时要注意避免因变量持久性等问题导致的错误，以确保数据分析的准确性和可靠性。

以下是广义最小二乘法步骤的mermaid流程图：

graph LR
    A[用文字描述问题] --> B[整理成标准形式]
    B --> C[检查数据]
    C --> D[确定数据准确性]
    D --> E[陈述先验信息]
    E --> F[整理先验信息]
    F --> G[确定先验信息准确性]
    G --> H[估计模型参数及协方差]
    H --> I[陈述估计值及置信区间]
    I --> J[检查个体误差]
    J --> K[检查总误差]
    K --> L[评估零假设（总误差）]
    M[有两个不同模型?] --> N[评估零假设（两个模型E值）]

以下是计算功率谱密度步骤的mermaid流程图：

graph LR
    A[计算时间参数] --> B[计算频率参数]
    B --> C[预处理时间序列]
    C --> D[计算傅里叶变换]
    D --> E[保留非负频率]
    E --> F[计算功率谱密度]
    F --> G[检查功率谱密度]
    G --> H[陈述零假设]
    H --> I[计算95%置信水平]
    I --> J[绘制置信水平并评估峰值显著性]

数据处理与分析中的关键方法及注意事项

7. 问题实践

在实际应用中，我们会遇到各种具体问题，下面列举了一些常见问题及解决步骤。

7.1 比较不同年份最热天均值

以黑岩森林温度数据集为例，该数据集的第一年和第十二年数据较为完整。我们要计算这两年中10个最热天（以最高峰值温度定义）的均值，并检验这两个均值是否有显著差异，步骤如下：
1. 陈述零假设 ：假设这两年10个最热天的均值没有显著差异。
2. 计算t统计量 ：根据数据计算该情况下的t统计量。
3. 判断是否拒绝零假设 ：根据t统计量的值判断是否能拒绝零假设。

7.2 分析预测误差滤波器的误差减少显著性

回顾Neuse河预测误差滤波器的计算，分析不同长度滤波器对之间误差减少的显著性。

7.3 证明变量的均值和方差特性

证明变量Z = (d - d’) / σd 的均值为零，方差为1，假设d服从均值为d’，方差为σd的正态分布，通过将概率密度函数p(d)转换为p(Z)来证明。

7.4 分析随机时间序列功率谱密度的峰值

对于随机时间序列的功率谱密度图（如图所示），完成以下操作：
1. 统计显著峰值数量 ：统计达到或超过95%显著水平的峰值数量，并与预期数量进行比较。
2. 计算最高峰值的显著水平 ：已知该数据集N = 1024，c = 0.634，计算最高峰值的显著水平。

7.5 确定Neuse河水位图年度峰值的显著性

确定Neuse河水位图功率谱密度中年度峰值的显著性，可参考相关脚本（如edama_06_13和edapy_06_13）。

7.6 检测叠加正弦波的置信度

修改脚本（如edama_12_13或edapy_12_13），在时间序列y上叠加小振幅A和频率f0的正弦波，评估在选定A和f0值下检测该正弦波的置信度。确定在检测置信度变差之前A可以小到什么程度，以y方差的平方根的分数形式表示。并探讨对于固定的A，高频或低频正弦波哪个更容易检测。

7.7 自助法统计斜率与截距的比率

修改脚本（如edama_12_14或edapy_12_14），生成斜率b与截距a的比率r = b / a的自助法统计数据。使用最小二乘法估计截距和斜率及其方差，并结合相关方法（如第11.5节的方法）生成可与自助法结果进行比较的解析结果。

8. 操作总结

为了更清晰地展示各种操作的关键信息，下面以表格形式总结了广义最小二乘法和计算功率谱密度的主要步骤。

操作类型	主要步骤
广义最小二乘法	1. 用文字描述问题；2. 整理成标准形式；3. 检查数据；4. 确定数据准确性；5. 陈述先验信息；6. 整理先验信息；7. 确定先验信息准确性；8. 估计模型参数及协方差；9. 陈述估计值及置信区间；10. 检查个体误差；11. 检查总误差；12. 若有两个模型，评估零假设（两个模型E值）
计算功率谱密度	1. 计算时间参数；2. 计算频率参数；3. 预处理时间序列；4. 计算傅里叶变换；5. 保留非负频率；6. 计算功率谱密度；7. 检查功率谱密度；8. 陈述零假设；9. 计算95%置信水平；10. 绘制置信水平并评估峰值显著性

同时，我们可以用以下mermaid流程图来直观展示问题实践的流程：

graph LR
    A[选择问题类型] --> B{不同年份最热天均值比较}
    A --> C{预测误差滤波器误差减少显著性分析}
    A --> D{证明变量均值和方差特性}
    A --> E{分析功率谱密度峰值}
    A --> F{确定年度峰值显著性}
    A --> G{检测叠加正弦波置信度}
    A --> H{自助法统计比率}
    B --> B1[陈述零假设]
    B1 --> B2[计算t统计量]
    B2 --> B3[判断是否拒绝零假设]
    C --> C1[分析滤波器误差减少情况]
    D --> D1[转换概率密度函数]
    E --> E1[统计显著峰值数量]
    E1 --> E2[计算最高峰值显著水平]
    F --> F1[参考相关脚本分析]
    G --> G1[修改脚本叠加正弦波]
    G1 --> G2[评估检测置信度]
    H --> H1[修改脚本生成统计数据]
    H1 --> H2[生成解析结果比较]

9. 总结与建议

在数据处理与分析过程中，我们介绍了多种关键方法，包括功率谱密度与AR1过程匹配、自助法置信区间估计、广义最小二乘法、变量管理和时间格式转换等。这些方法在不同的数据分析场景中都有重要应用，但也需要注意一些问题。

• 变量管理 ：在MATLAB和Python中，变量的持久性可能会导致一些错误，如使用错误的变量值、脚本意外失败等。因此，在运行脚本前，建议删除工作区中的所有变量，避免因变量残留导致的问题。
• 时间格式转换 ：对于包含日历和时钟格式时间的数据集，要及时将其转换为适合数据分析的格式。根据不同的分析需求，选择合适的时间单位和起始时间。
• 方法选择 ：在处理复杂的数据分析问题时，如确定具有复杂关系的模型参数的置信区间，自助法是一种有效的方法。而广义最小二乘法适用于解决数据与模型关系的估计问题。

通过合理运用这些方法，并注意相关问题，我们可以提高数据分析的准确性和可靠性，更好地从数据中提取有价值的信息。

希望以上内容能帮助你在数据处理与分析中更加得心应手，如果你在实践过程中遇到问题，可以随时参考上述方法和步骤进行解决。