计算物理（八）数理统计
fengxiaot Lv4
  2023-05-05 21:36:18 2023-05-05 21:36  

本文介绍如何用 Python 实现文件I/O、插值、拟合、随机数。

文件I/O

Python 提供了必要的函数和方法进行默认情况下的文件基本操作。你可以用 file 对象做大部分的文件操作。

打开文件

Python 内置 open() 函数来打开文件，并创建一个 file 对象

1
file = open(path, mode="r")
  • path 为字符串，包含文件的路径和名称
  • mode 指定文件的打开模式
模式 描述
r 以只读方式打开文件。文件指针将会放在文件的开头。
r+ 打开一个文件用于读写。文件指针将会放在文件的开头，原有内容不会被预先清空 (truncate) 。
w 打开一个文件并清空所有内容用于写入。文件指针将会放在文件的开头。如果该文件不存在，创建新文件。
w+ 打开一个文件用于读写。文件指针将会放在文件的开头，原有内容会被预先清空。如果该文件不存在，创建新文件。
a 打开一个文件用于追加写。文件指针将会放在文件的结尾。如果该文件不存在，创建新文件进行写入。
a+ 打开一个文件用于追加读写。文件指针将会放在文件的结尾。如果该文件不存在，创建新文件用于读写。

关闭文件

调用 file 对象的 close() 方法以关闭文件。

Python 为防止程序员忘记关闭文件，提供 with 关键字以打开文件，并在语句执行完毕后自动关闭文件。

1
2
3
4
5
6
# Method 1
file = open(path, mode="r")
file.close()
# Method 2
with open(path, mode="r") as file:
    # Operations to the file

with 关键字实际上等价于 try/finally 语句，with 关键字支持用 as 创建别名。

文件写入

Python 文件写入操作使用 file 对象的 write() 方法或 writelines() 方法

write()

  • write() 方法用于向文件中写入单个字符串。
  • 它接受一个字符串作为参数，并将该字符串写入文件。
  • 如果多次调用 write() 方法，每次调用将在文件中追加写入的内容。

writelines()

  • writelines() 方法用于向文件中写入多个字符串。
  • 它接受一个可迭代对象（如列表）作为参数，其中每个元素都是字符串。
  • 它将迭代对象中的每个字符串写入文件。

/* 注：无论 write() 还是 writelines() 都不会自动向字符串末尾添加换行符’\n’，需要手动插入，不要被 writelines() 中的 line 所迷惑！ */

文件读取

Python 文件读取使用 file 对象的 read() 方法、readline() 方法或 readlines() 方法

1
file.read(size)

其中 size 为读取的字节数，默认为 -1 ，表示读取全部内容。转义字符 '\n' 等也占且仅占 1 个字节。

假设 foo.txt 包含以下内容

1
2
Strawberry Creek Park
Berkeley, CA, 94704

执行如下 Python 代码

1
2
3
str1 = file.read(10)       # str1 = 'Strawberry'
str2 = file.readline()     # str2 = 'Strawberry Creek Park\n'
strlist = file.readlines() # str3 = ['Strawberry Creek Park\n','Berkeley, CA, 94704']

注：read(), readline() 和 readlines() 读取的字符串都包含行末换行符！

可以使用 .rstrip() 字符串方法删除字符串末尾的空白字符，如空格、制表符、换行符等。

文件操作

清空内容 使用 file 对象的 truncate() 方法

光标定位 使用 file 对象的 tell() 方法获知光标的当前位置，使用 seek(offset, whence=0) 移动光标


CSV 数据

对 CSV 数据的操作可使用 Pandas ，一个强大的分析结构化数据的工具集

1
import pandas as pd

DataFrame

DataFrame 是 Pandas 内建的表格型数据结构，包含一组带标签的列，每列可以是不同的值类型（数值、字符串、布尔型值），并共用一个行索引，默认为 RangeIndex(0,…,n) 。

1
class pandas.DataFrame(data, index=None, columns=None, dtype=None, copy=None)
  • data 为数据，可从 dict 或 ndarray 构建
  • index 为行索引，默认为 0,1,⋯ ,n−10,1,\cdots,n-1
  • columns 为列索引，可指定为 list of strings 标签，方便说明每列数据的意义。默认为 0,1,⋯ ,n−10,1,\cdots,n-1 。
  • dtype 为数据类型，默认为 None 自动推断
  • copy 指定是否复制一份源数据的副本，若为 None 则依据内置规则确定

如下数据可创建为 DataFrame 对象

stu age sex
0 ‘Feng’ 20 ‘M’
1 ‘Yaki’ 21 ‘F’
2 ‘J-20’ 7 ‘J’ 
1
2
3
4
5
6
# Method 1
data = {'stu': ['Feng','Yaki','J-20'], 'age': [20,21,7], 'sex':['M','F','J']}
df = pandas.DataFrame(data)
# Method 2
list2d = [['Feng',20,'M'],['Yaki',21,'F'],['J-20',7,'J']]
df = pandas.DataFrame(list2d, columns=['stu','age','sex'])

文件读取

外置 CSV 文件的读取使用 read_csv 函数

1
pandas.read_csv(file, header=None, names=no_default)
  • header 指示是否包含列标签，默认为 infer
  • names 指定列标签，为字符串列表

数据索引

列数据 使用 dfobj['col'].to_numpy() 获取列标签为 'col' 的转为 numpy 数据类型后的数据


插值

插值使用 SciPy 中的 interpolate 模块。

一维插值

插值用于已知若干离散数据点 (xi,yi)(x_i, y_i)，在其间构造连续函数 f(x)f(x)。

1
2
f = interp1d(x, y, kind='linear')
y0 = f(x0)

interp1d 的 kind 参数可指定插值方法，默认为线性插值 (linear) ，还可指定为最近邻插值 (nearest) 、二次样条插值 (quadratic)，三次样条插值 (cubic) 。一般使用 cubic 。

interp1d 返回一个可调用的函数，而非离散的预测值。

interp1d 不支持外推，若传入 x_new 超过原始区间，会报错。

多维插值

二维或三维插值通常使用以下函数：

  • scipy.interpolate.griddata：任意点分布的插值，将散点转化为网格
  • scipy.interpolate.RectBivariateSpline：规则网格的二维样条插值
1
2
3
4
5
6
7
8
from scipy.interpolate import griddata

points = np.random.rand(100, 2)
values = np.sin(points[:,0]) + np.cos(points[:,1])

grid_x, grid_y = np.mgrid[0:1:100j, 0:1:100j]

grid_z = griddata(points, values, (grid_x, grid_y), method='cubic')

拟合

拟合使用 SciPy 中的 Optimize 模块。

拟合的典型任务为给定观测数据 (xi,yi)(x_i, y_i)，求模型函数 f(x,θ)f(x,\theta) 的最优参数 θ\theta。

一维非线性拟合

一维非线性拟合使用 curve_fit

1
from scipy.optimize import curve_fit

基本语法

1
popt, pcov = curve_fit(func, xdata, ydata, p0=None)

其中

  • func：模型函数，如 f(x,a,b)=a*x+b
  • popt：最优拟合参数
  • pcov：参数协方差矩阵

示例：拟合 y=Ae−kxy = A e^{-kx} 形式的衰减函数

1
2
3
4
5
6
7
8
9
10
11
import numpy as np
from scipy.optimize import curve_fit

f = lambda x, A, k: A * np.exp(-k * x)

x = np.linspace(0, 5, 50)
y = 3 * np.exp(-1.2 * x) + 0.1*np.random.randn(50)

popt, pcov = curve_fit(f, x, y)

A_fit, k_fit = popt

多维拟合与一般优化

更复杂的目标函数可使用 scipy.optimize.minimize

1
2
cost = lambda p: np.sum((y - f(x, p))**2)
res = minimize(cost, p0)

随机数

NumPy 提供 random 模块用于生成各种分布的随机数。

1
rng = np.random.default_rng(seed=42)

基本分布

1
2
3
4
5
6
rng.random(size)              # 均匀分布 U(0,1)
rng.normal(mu, sigma, size)   # 正态分布 N(mu, sigma)
rng.integers(low, high, size) # 整数随机数
rng.poisson(lam, size)        # 泊松分布
rng.exponential(scale, size)  # 指数分布
rng.uniform(a, b, size)       # 区间均匀分布

多维随机矩阵

1
2
rng.random((3,3))     # 3x3 矩阵
rng.normal(0,1,(5,5)) # 高斯随机矩阵