超级会员免费看
环境数据分析中的Python与MATLAB实用技巧
1. 时间算术与闰秒问题
在进行时间算术运算时，Python 可以轻松计算时间间隔。例如：
# edapy_13_02: time arithmetic
from datetime import datetime
DELTA = datetime(2008, 2, 11, 4, 4, 1) - datetime(2008, 2, 11, 4, 4, 0)
print('DT: %.2f\n' % DELTA.total_seconds())
不过，当需要精确到秒或更高精度的时间计算时，闰秒是一个复杂的问题。闰秒类似于闰年，是对时钟的整秒修正，通常在每年的6月30日和12月31日进行，用于补偿地球自转的微小不规则性。与闰年不同，闰秒由国际地球自转和参考系统服务局（IERS）每半年确定一次。
这就导致跨越6月30日和12月31日的时间间隔，如果没有最新的闰秒表，就无法准确计算。更糟糕的是，广泛使用的协调世界时（UTC）使用闰秒，而像全球定位系统（GPS）等其他标准则不使用。截至2020年底，自1972年首次实施闰秒以来，总共宣布了34个闰秒。如果不考虑闰秒，长达数十年的时间间隔可能会产生数十秒的误差。需要注意的是，MATLAB的
datenum()
函数和Python的
datetime()
方法都没有考虑闰秒，因此对于UTC时间无法提供秒级精度。
2. 读取复杂文本文件
MATLAB的
load()
函数和Python的
np.genfromtxt()
方法只能读取包含数值表的文本文件。但一些公开数据库提供的是更复杂的文本文件，包含数字和字母混合的值。例如，Black Rock Forest温度数据集包含时间和温度信息，文件内容如下：
2100-2159 31 Jan 1997
-1.34
2200-2259 31 Jan 1997
-0.958
2300-2400 31 Jan 1997
-0.601
0000-0059 1 Feb 1997
-0.245
0100-0159 1 Feb 1997
-0.217
对于这种简单格式的文件，可以使用Microsoft Excel电子表格软件的文本导入向导模块进行重新格式化。但对于更复杂的情况，目前没有通用且易于使用的软件能够可靠处理，通常需要为每个文件编写自定义脚本。
以下是MATLAB和Python读取此类文件的示例代码：
MATLAB代码
：
% edapy_13_03: reading idiosyncratic text file
filename = '../Data/brf_raw.txt';
fid = fopen(filename, "r");
N = 0; % line number counter
while 1 % loop over lines in file
myline = fgetl(fid); % read line
if ~ischar(myline)
% check for end of file
fclose(fid);
break;
end
N = N + 1; % increment line counter
% now process the line
% ...
end
Python代码 ：
# edapy_13_03: reading idiosynchratic text file
filename = '../Data/brf_raw.txt';
fid = open(filename, "r");
N = 0; # line number counter
while 1: # loop over lines in file
myline = fid.readline(); # read line
if not myline:
fid.close();
break;
# now process the line
# ...
N = N + 1;
在处理文件时，首先需要打开文件，使用变量
fid
引用打开的文件。然后在循环中逐行读取文件内容到字符串
myline
中，当读取到最后一行时循环终止。变量
N
用于统计行数。接下来需要对字符串进行处理，将所有数据字段转换为数值并存储在数组中。最后，关闭文件。
在处理部分，MATLAB的
replace()
、
strsplit()
和
str2num()
函数非常有用，Python中对应的方法有
myline.replace()
、
myline.split()
、
int()
和
float()
。数据格式转换脚本编写起来很繁琐，一定要仔细测试，包括与原始数据进行抽查比对，抽查应包括文件末尾附近的数据。
3. 误差传播规则
假设使用线性规则从N个数据
d
中形成MA个模型参数
m(A)
，即
m(A) = M(A)d
。当
MA = N
时，协方差矩阵的关系为
Cm(A) = MCdMT
。为了验证
MA < N
的情况，首先设计
MB = N - MA
个互补模型参数
m(B)
，使得
m(B) = M(B)d
。然后将两组模型参数连接起来，得到一个方阵方程：
[
\begin{bmatrix}
m(A) \
m(B)
\end{bmatrix}
=
\begin{bmatrix}
M(A) \
M(B)
\end{bmatrix}
d
]
根据正常的误差传播规则，可得：
[
C_m =
\begin{bmatrix}
M(A) \
M(B)
\end{bmatrix}
C_d
\begin{bmatrix}
M(A)^T & M(B)^T
\end{bmatrix}
=
\begin{bmatrix}
M(A)C_dM(A)^T & M(A)C_dM(B)^T \
M(B)C_dM(A)^T & M(B)C_dM(B)^T
\end{bmatrix}
=
\begin{bmatrix}
C_m^{(A)} & C_m^{(A,B)} \
C_m^{(B,A)} & C_m^{(B)}
\end{bmatrix}
]
其中，
Cm(A)
包含了
m(A)
模型参数的所有方差和协方差，满足正常的误差传播规则，并且与
M(B)
的选择无关。因此，该规则可以应用于
MA < N
且
M(A)
为矩形的情况，无需考虑互补模型参数的具体选择。
4. eda_draw()函数
可以使用
eda_draw()
函数将一系列方阵和向量绘制成灰度图像，还可以在矩阵和向量下方添加标题，并在它们之间绘制符号。以下是MATLAB和Python的示例代码：
MATLAB代码
：
% edama_13_04, example of use of eda_draw() function
% ...
eda_draw(d, 'caption d', '=', G, 'caption G', m, 'caption m');
Python代码 ：
# edapy_13_04, example of use of eda_draw() function
# ...
eda_draw( 'title d', d, '=', 'title G', G, 'title m', m);
该函数可以创建方程
d = Gm
的图形表示，它可以接受向量、方阵和字符串，字符串可以放在图像之间和下方。
5. 复数最小二乘法
在处理包含复数的最小二乘法问题时（例如，当模型参数是函数的傅里叶变换时），方程
Gm = d
中的所有量都是复数。总误差
E
的正确定义为：
[
E(m) = e^T
e = (d - Gm)^
T(d - Gm)
]
其中
*
表示复共轭，这种复共轭和矩阵转置的组合称为埃尔米特转置，记为
e^H = e^T*
。总误差
E
是一个非负实数。通过将
m
表示为
m = m(R) + im(I)
，并对实部和虚部分别求偏导数并令其为零，可以得到最小二乘解和协方差：
[
m_{est} = (G^HG)^{-1}G^Hd
]
[
C_m = \sigma_d^2(G^HG)^{-1}
]
在MATLAB中，复数矩阵
G
的埃尔米特转置与转置使用相同的符号，如
G'
，无复共轭的转置为
G.'
，因此实现复数最小二乘法时无需更改MATLAB公式。在Python中，矩阵
G
的埃尔米特转置计算为
G.H
。
6. 广义最小二乘法的推导
严格来说，在方程
Hm = h
中，只有当
K = M
（即矩阵
H
为方阵且
H^-1
存在）时，长度为
K
的
h
的概率密度函数
p(h)
才能与长度为
M
的
m
的概率密度函数
p(m)
成比例。在其他情况下，雅可比行列式是未定义的。
当只有少量先验信息可用时，会出现非方阵的情况。可以通过添加
M - K
行互补信息使
H
成为方阵，并赋予它们可忽略的确定性，从而不影响广义最小二乘解。这种修正不会影响推导结果，广义最小二乘解及其协方差的所有公式保持不变。根本问题在于，代表无信息状态的均匀概率密度函数在无界域上不存在，最好的方法是使用非常宽的正态概率密度函数。
7. MATLAB和Python函数
MATLAB和Python都可以定义类似于内置函数（如
sin()
和
cos()
）的自定义函数。以下是定义计算圆面积的函数示例：
MATLAB代码
：
function a = areaofcircle(r)
% computes area, a, of circle of radius, r.
a = pi * (r^2);
return
end
Python代码 ：
# edapy_13_05, simple examples of functions
# computes area, a, of circle of radius, r.
def areaofcircle(r):
a = pi*(radius**2);
return a;
虽然具体语法不同，但MATLAB和Python函数有很多共同特点。函数都有名称（如
areaofcircle
），输入参数（如
r
）在名称后的括号中指定，用于函数内部的计算。计算结果（如
a
）在函数调用时返回。
在MATLAB中，函数定义必须放在LiveScript的最后一个单元格中；在Python中，可以放在任何位置，通常放在Jupyter Notebook的第一个单元格中。调用函数的示例如下：
MATLAB代码
：
% edama_13_05: simple examples of a function
radius = 2;
area = areaofcircle(radius);
Python代码 ：
# edapy_13_05, simple examples of functions
radius = 2;
area = areaofcircle(radius);
需要注意的是，主脚本中的变量名不必与函数中的变量名匹配，函数中的变量只是占位符。
函数还可以接受多个输入变量并返回多个输出变量，例如计算矩形周长和面积的函数：
MATLAB代码
：
function [c,a] = CandAofrectangle(l, w)
% computes circumference, c, and area, a, of
% a rectangle of length, l, and width, w.
c = 2*(l+w);
a = l*w;
return
end
Python代码 ：
# edapy_13_05, simple examples of functions
# computes circumference, c, and area, a, of
# a rectangle of length, l, and width, w
def CandAofRectangle(a,b):
circ = 2*(a+b);
area = a*b;
return circ, area;
调用示例：
MATLAB代码
：
% edama_13_05: simple examples of a function
a = 2;
b = 4;
[circ, area] = CandAofrectangle(a,b);
Python代码 ：
# edapy_13_05, simple examples of functions
a = 2;
b = 4;
circ, area = CandAofRectangle(a,b);
8. 矩阵重组
在MATLAB和Python中，只要
N1M1 = N2M2
，就可以将
N1 x M1
的矩阵形状更改为
N2 x M2
的矩阵。例如，一个
4 x 4
的矩阵可以轻松转换为等效的
1 x 16
、
2 x 8
、
8 x 2
和
16 x 1
矩阵。以下是将
3 x 4
矩阵重组为
12 x 1
列向量的示例：
MATLAB代码
：
% eda13_06: examples of reorganizing a matrix
A = [ [1,2,3,4]; [5,6,7,8]; [9,10,11,12] ];
[N,M] = size(A);
% create a column-vector, c, containing the elements of A,
% organized column-wise.
c = reshape(A,M*N,1);
Python代码 ：
# edapy_13_06, examples of reorganizing a matrix
import numpy as np
A = np.array( [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]] );
N,M = np.shape(A);
# create a column-vector, b, containing the elements of A,
# organized column-wise.
c = np.reshape(A,(N*M,1),order='F');
这些函数通常可以消除脚本中矩阵重组部分的循环。在MATLAB中，矩阵元素可以使用单个索引访问，而不是通常的两个索引。例如，矩阵
A
通常用
A(i,j)
寻址，也可以用
A(k)
寻址，其中
k
是按列展开
A
的所有元素形成的列向量的索引。这种寻址方式可以简化某些操作，但建议在标准寻址过于繁琐时谨慎使用。
9. atan2()函数的使用
离散傅里叶变换的相位
φ = tan^-1(B/A)
定义在区间
(-π, +π)
上。在脚本中，应该使用
atan2(B,A)
函数，而不是
atan(B/A)
函数。后者定义在错误的区间
(-π/2, +π/2)
上，并且当
A = 0
时会失败。
10. 离散傅里叶数据核的正交规范性
对于复数版本的傅里叶数据核
G
，规则
[G^T*G] = NI
可以推导如下。首先，写出傅里叶级数未归一化的数据核定义：
[
G_{kp} = exp(2\pi i(k - 1)(p - 1)/N)
]
则
[G^T*G]
为：
[
(G^T
G)
{pq} = \sum
{k = 1}^{N} G_{kp}^
G_{kq} = \sum_{k = 1}^{N} exp(2\pi i k(p - q)/N) = \sum_{k = 0}^{N - 1} z^k = f(z)
]
其中
f(z) = exp(2\pi i(p - q)/N)
。考虑级数
f(z) = \sum_{k = 0}^{N - 1} z^k = 1 + z + z^2 + \cdots + z^{N - 1}
，乘以
z
得到
zf(z) = \sum_{k = 0}^{N - 1} z^{k + 1} = z + z^2 + \cdots + z^N
，相减可得：
[
f(z) - zf(z) = 1 - z^N
]
即：
[
f(z) = \frac{1 - z^N}{1 - z}
]
将
z = exp(2\pi i(p - q)/N)
代入，当
p ≠ q
时，分子为零，分母不为零，所以
f(z) = 0
，即
[G^T*G]
的非对角元素为零。当
p = q
时，分母也为零，使用洛必达法则求极限可得：
[
f(z) = \lim_{s \to 0} \frac{2\pi i exp(2\pi i s)}{(1/N) exp(2\pi i s/N)} = N
]
所以
[G^T*G]
的对角元素都等于
N
。
11. 函数在正交基下的展开
假设在区间
(tmin, tmax)
上用一组基函数
gi(t)
的和来近似任意函数
d(t)
：
[
d(t) \approx d(s)
]
[
d(s) = \sum_{i = 1}^{M} m_i g_i(t)
]
其中
mi
是未知系数。傅里叶级数就是这样一种近似，基函数为正弦和余弦函数：
[
g_i(t) =
\begin{cases}
\cos(\omega_j t) & i \text{ 为奇数} \
\sin(\omega_j t) & i \text{ 为偶数}
\end{cases}
]
其中
j = floor((i + 1)/2)
。对于任何给定的系数集，可以通过定义误差来衡量近似的质量：
[
E = \int_{t_{min}}^{t_{max}} (d(t) - d(s))^2 dt
]
这是通常最小二乘误差的推广，当
d(s) = d(t)
时，
E = 0
。一般来说，只有当
M → ∞
时才能实现零误差。通过将
E
视为未知系数的函数，并对其求偏导数
∂E/∂mk = 0
，可以得到未知系数的方程：
[
\int_{t_{min}}^{t_{max}} g_j(t) d(t) dt = \sum_{i = 1}^{M} m_i \int_{t_{min}}^{t_{max}} g_j(t) g_i(t) dt
]
即
b = Mm
，其中：
[
b_j = \int_{t_{min}}^{t_{max}} g_j(t) d(t) dt
]
[
M_{ij} = \int_{t_{min}}^{t_{max}} g_j(t) g_i(t) dt
]
求解系数可得
m = M^-1b
。在某些情况下，基函数
gi(t)
可能是正交归一的，即：
[
\int_{t_{min}}^{t_{max}} g_j(t) g_i(t) dt = \delta_{ij}
]
其中
δij
是克罗内克δ函数。在这种情况下，
M = M^-1 = I
，系数公式简化为
m = b
，每个系数独立确定，与求和项的数量甚至其他基函数的身份无关。
环境数据分析中的Python与MATLAB实用技巧（续）
11. 函数在正交基下的展开（续）
在前面我们提到了函数在正交基下展开的相关内容，这里我们进一步梳理其流程，方便大家理解和操作。以下是函数在正交基下展开的步骤：
1.
定义近似函数
：假设要在区间
(tmin, tmax)
上用一组基函数
gi(t)
的和来近似任意函数
d(t)
，表示为：
[
d(t) \approx d(s)
]
[
d(s) = \sum_{i = 1}^{M} m_i g_i(t)
]
这里的
mi
是未知系数。
2.
选择基函数
：例如傅里叶级数的基函数为正弦和余弦函数：
[
g_i(t) =
\begin{cases}
\cos(\omega_j t) & i \text{ 为奇数} \
\sin(\omega_j t) & i \text{ 为偶数}
\end{cases}
]
其中
j = floor((i + 1)/2)
。
3.
定义误差函数
：通过定义误差来衡量近似的质量，误差函数为：
[
E = \int_{t_{min}}^{t_{max}} (d(t) - d(s))^2 dt
]
当
d(s) = d(t)
时，
E = 0
。一般只有当
M → ∞
时才能实现零误差。
4.
求未知系数方程
：将
E
视为未知系数的函数，对其求偏导数
∂E/∂mk = 0
，得到未知系数的方程：
[
\int_{t_{min}}^{t_{max}} g_j(t) d(t) dt = \sum_{i = 1}^{M} m_i \int_{t_{min}}^{t_{max}} g_j(t) g_i(t) dt
]
即
b = Mm
，其中：
[
b_j = \int_{t_{min}}^{t_{max}} g_j(t) d(t) dt
]
[
M_{ij} = \int_{t_{min}}^{t_{max}} g_j(t) g_i(t) dt
]
5.
求解系数
：求解上述方程得到
m = M^-1b
。如果基函数是正交归一的，即：
[
\int_{t_{min}}^{t_{max}} g_j(t) g_i(t) dt = \delta_{ij}
]
其中
δij
是克罗内克δ函数，此时
M = M^-1 = I
，系数公式简化为
m = b
。
下面用一个简单的mermaid流程图来表示这个过程：
graph TD;
A[定义近似函数d(t)≈d(s)] --> B[选择基函数gi(t)];
B --> C[定义误差函数E];
C --> D[求偏导数∂E/∂mk = 0];
D --> E[得到方程b = Mm];
E --> F{基函数是否正交归一};
F -- 是 --> G[m = b];
F -- 否 --> H[m = M^-1b];
总结
在环境数据分析中，Python和MATLAB提供了丰富的工具和函数来处理各种问题。我们从时间算术运算开始，了解了闰秒对高精度时间计算的影响，以及MATLAB和Python在处理时间相关问题时的局限性。在读取复杂文本文件方面，虽然有简单情况可以借助Excel的文本导入向导，但复杂情况仍需编写自定义脚本，同时介绍了处理脚本中有用的函数和方法。误差传播规则的推导让我们明白在不同情况下如何处理模型参数的协方差。
eda_draw()
函数为矩阵和向量的可视化提供了便利。复数最小二乘法和广义最小二乘法的推导则为处理复杂的数据分析问题提供了理论支持。
在函数定义和使用方面，MATLAB和Python都允许用户自定义函数，并且函数可以接受多个输入和返回多个输出。矩阵重组功能让我们可以灵活改变矩阵的形状，提高代码的效率。
atan2()
函数的正确使用确保了离散傅里叶变换相位计算的准确性。离散傅里叶数据核的正交规范性推导加深了我们对傅里叶变换的理解。最后，函数在正交基下的展开为函数的近似表示提供了一种有效的方法。
以下是一个总结表格，对比了MATLAB和Python在不同方面的使用：
| 功能 | MATLAB | Python |
| ---- | ---- | ---- |
| 时间算术 |
datenum()
不考虑闰秒 |
datetime()
不考虑闰秒 |
| 读取复杂文本文件 |
fopen()
、
fgetl()
、
replace()
、
strsplit()
、
str2num()
|
open()
、
readline()
、
myline.replace()
、
myline.split()
、
int()
、
float()
|
| 误差传播规则 | 按规则推导计算 | 按规则推导计算 |
| 绘图函数 |
eda_draw()
|
eda_draw()
|
| 复数最小二乘法 | 用
G'
计算埃尔米特转置 | 用
G.H
计算埃尔米特转置 |
| 函数定义 | 定义在LiveScript最后单元格 | 可定义在Jupyter Notebook任意位置，常放第一个单元格 |
| 矩阵重组 |
reshape()
，可单索引访问元素 |
np.reshape()
|
| 相位计算 | 用
atan2()
| 用
atan2()
|
通过掌握这些实用技巧，我们可以更高效地进行环境数据分析，解决实际问题。在实际应用中，我们可以根据具体需求选择合适的工具和方法，同时要注意代码的测试和验证，确保结果的准确性。希望这些内容对大家在环境数据分析的学习和实践中有所帮助。
扫一扫
427
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
