SciPy 基础概念

在上一章中，我们学习了如何安装和配置 SciPy。本章将深入探讨 SciPy 的基础概念，包括数据结构、常用函数和编程模式。

SciPy 的核心理念

1. 基于 NumPy 数组

SciPy 的所有功能都建立在 NumPy 数组之上。理解 NumPy 数组是使用 SciPy 的基础：

python

import numpy as np
import scipy as sp

# NumPy 数组是 SciPy 的基础
data = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
print(f"数据类型: {type(data)}")
print(f"数组形状: {data.shape}")
print(f"数据类型: {data.dtype}")

2. 模块化设计

SciPy 采用模块化设计，每个模块专注于特定的科学计算领域：

python

# 导入特定模块
from scipy import stats, optimize, integrate, linalg

# 或者导入整个 scipy 包
import scipy

# 查看可用的子模块
print(dir(scipy))

3. 函数式编程风格

SciPy 主要采用函数式编程风格，大多数操作通过函数调用完成：

python

import numpy as np
from scipy import stats

# 函数式调用
data = np.random.normal(0, 1, 100)
mean = np.mean(data)
std_dev = np.std(data)
t_stat, p_value = stats.ttest_1samp(data, 0)

print(f"均值: {mean:.3f}")
print(f"标准差: {std_dev:.3f}")
print(f"t 检验 p 值: {p_value:.3f}")

常用数据类型

1. 标量（Scalars）

python

import numpy as np

# 不同类型的标量
integer_scalar = np.int32(42)
float_scalar = np.float64(3.14159)
complex_scalar = np.complex128(1 + 2j)

print(f"整数: {integer_scalar}, 类型: {type(integer_scalar)}")
print(f"浮点数: {float_scalar}, 类型: {type(float_scalar)}")
print(f"复数: {complex_scalar}, 类型: {type(complex_scalar)}")

2. 一维数组（Vectors）

python

# 创建一维数组
vector = np.array([1, 2, 3, 4, 5])
print(f"向量: {vector}")
print(f"形状: {vector.shape}")
print(f"维度: {vector.ndim}")

# 常用的一维数组创建方法
zeros_vec = np.zeros(5)
ones_vec = np.ones(5)
range_vec = np.arange(0, 10, 2)
linspace_vec = np.linspace(0, 1, 5)

print(f"零向量: {zeros_vec}")
print(f"单位向量: {ones_vec}")
print(f"等差数列: {range_vec}")
print(f"等分数列: {linspace_vec}")

3. 二维数组（Matrices）

python

# 创建二维数组
matrix = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(f"矩阵:\n{matrix}")
print(f"形状: {matrix.shape}")
print(f"维度: {matrix.ndim}")

# 常用的二维数组创建方法
zeros_mat = np.zeros((3, 3))
ones_mat = np.ones((2, 4))
identity_mat = np.eye(3)
random_mat = np.random.random((2, 3))

print(f"零矩阵:\n{zeros_mat}")
print(f"单位矩阵:\n{identity_mat}")
print(f"随机矩阵:\n{random_mat}")

4. 多维数组（Tensors）

python

# 创建三维数组
tensor = np.random.random((2, 3, 4))
print(f"张量形状: {tensor.shape}")
print(f"张量维度: {tensor.ndim}")
print(f"张量大小: {tensor.size}")

# 访问张量元素
print(f"第一个 2D 切片:\n{tensor[0]}")
print(f"特定元素: {tensor[0, 1, 2]}")

基本数组操作

1. 数组索引和切片

python

import numpy as np

# 一维数组索引
arr_1d = np.array([10, 20, 30, 40, 50])
print(f"第一个元素: {arr_1d[0]}")
print(f"最后一个元素: {arr_1d[-1]}")
print(f"切片 [1:4]: {arr_1d[1:4]}")

# 二维数组索引
arr_2d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(f"元素 [1,2]: {arr_2d[1, 2]}")
print(f"第二行: {arr_2d[1, :]}")
print(f"第三列: {arr_2d[:, 2]}")

# 布尔索引
data = np.array([1, 2, 3, 4, 5, 6])
mask = data > 3
print(f"大于 3 的元素: {data[mask]}")

2. 数组形状操作

python

# 改变数组形状
arr = np.arange(12)
print(f"原始数组: {arr}")

# 重塑为 2D 数组
reshaped = arr.reshape(3, 4)
print(f"重塑后:\n{reshaped}")

# 展平数组
flattened = reshaped.flatten()
print(f"展平后: {flattened}")

# 转置
transposed = reshaped.T
print(f"转置后:\n{transposed}")

3. 数组拼接和分割

python

# 数组拼接
arr1 = np.array([1, 2, 3])
arr2 = np.array([4, 5, 6])

# 水平拼接
hstacked = np.hstack([arr1, arr2])
print(f"水平拼接: {hstacked}")

# 垂直拼接
vstacked = np.vstack([arr1, arr2])
print(f"垂直拼接:\n{vstacked}")

# 数组分割
data = np.arange(10)
split_arrays = np.split(data, 5)
print(f"分割结果: {split_arrays}")

SciPy 常用函数模式

1. 统计函数模式

python

from scipy import stats
import numpy as np

# 生成测试数据
np.random.seed(42)
data = np.random.normal(100, 15, 1000)

# 描述性统计
print(f"均值: {np.mean(data):.2f}")
print(f"中位数: {np.median(data):.2f}")
print(f"标准差: {np.std(data):.2f}")
print(f"偏度: {stats.skew(data):.3f}")
print(f"峰度: {stats.kurtosis(data):.3f}")

# 假设检验
t_stat, p_value = stats.ttest_1samp(data, 100)
print(f"单样本 t 检验: t={t_stat:.3f}, p={p_value:.3f}")

# 分布拟合
mu, sigma = stats.norm.fit(data)
print(f"拟合的正态分布参数: μ={mu:.2f}, σ={sigma:.2f}")

2. 优化函数模式

python

from scipy import optimize
import numpy as np

# 定义目标函数
def objective_function(x):
    return x**2 + 10*np.sin(x)

# 寻找最小值
result = optimize.minimize_scalar(objective_function)
print(f"最小值位置: x={result.x:.3f}")
print(f"最小值: f(x)={result.fun:.3f}")

# 求解方程
def equation(x):
    return x**3 - 2*x - 5

root = optimize.fsolve(equation, 2.0)[0]
print(f"方程的根: x={root:.3f}")
print(f"验证: f({root:.3f})={equation(root):.6f}")

3. 积分函数模式

python

from scipy import integrate
import numpy as np

# 定义被积函数
def integrand(x):
    return np.exp(-x**2)

# 数值积分
result, error = integrate.quad(integrand, 0, np.inf)
print(f"积分结果: {result:.6f}")
print(f"估计误差: {error:.2e}")
print(f"理论值: {np.sqrt(np.pi)/2:.6f}")

# 多重积分
def integrand_2d(y, x):
    return x * y

result_2d = integrate.dblquad(integrand_2d, 0, 1, lambda x: 0, lambda x: 1)
print(f"二重积分结果: {result_2d[0]:.6f}")

错误处理和调试

1. 常见错误类型

python

import numpy as np
from scipy import stats

# 维度不匹配错误
try:
    a = np.array([1, 2, 3])
    b = np.array([[1, 2], [3, 4]])
    result = a + b  # 这会引发广播错误
except ValueError as e:
    print(f"维度错误: {e}")

# 数值计算错误
try:
    result = np.log(-1)  # 负数的对数
except RuntimeWarning as e:
    print(f"数值警告: {e}")

# 统计函数错误
try:
    empty_data = np.array([])
    mean = np.mean(empty_data)
except RuntimeWarning as e:
    print(f"空数组警告: {e}")

2. 调试技巧

python

# 检查数组属性
def debug_array(arr, name="数组"):
    print(f"{name} 信息:")
    print(f"  形状: {arr.shape}")
    print(f"  数据类型: {arr.dtype}")
    print(f"  维度: {arr.ndim}")
    print(f"  大小: {arr.size}")
    print(f"  内存使用: {arr.nbytes} 字节")
    print(f"  是否连续: {arr.flags.c_contiguous}")
    print(f"  最小值: {np.min(arr)}")
    print(f"  最大值: {np.max(arr)}")
    print(f"  是否包含 NaN: {np.any(np.isnan(arr))}")
    print(f"  是否包含无穷大: {np.any(np.isinf(arr))}")
    print()

# 使用示例
test_array = np.random.random((3, 4))
debug_array(test_array, "测试数组")

性能优化基础

1. 向量化操作

python

import numpy as np
import time

# 避免使用循环，使用向量化操作
n = 1000000

# 低效的循环方式
start_time = time.time()
result_loop = []
for i in range(n):
    result_loop.append(i ** 2)
loop_time = time.time() - start_time

# 高效的向量化方式
start_time = time.time()
data = np.arange(n)
result_vectorized = data ** 2
vectorized_time = time.time() - start_time

print(f"循环方式耗时: {loop_time:.4f} 秒")
print(f"向量化耗时: {vectorized_time:.4f} 秒")
print(f"性能提升: {loop_time / vectorized_time:.1f} 倍")

2. 内存管理

python

# 使用合适的数据类型
data_float64 = np.random.random(1000000)  # 默认 float64
data_float32 = np.random.random(1000000).astype(np.float32)  # 转换为 float32

print(f"float64 内存使用: {data_float64.nbytes / 1024 / 1024:.2f} MB")
print(f"float32 内存使用: {data_float32.nbytes / 1024 / 1024:.2f} MB")
print(f"内存节省: {(1 - data_float32.nbytes / data_float64.nbytes) * 100:.1f}%")

# 就地操作
data = np.random.random(1000000)

# 创建新数组（消耗更多内存）
result1 = data * 2

# 就地操作（节省内存）
data *= 2  # 等价于 data = data * 2，但更节省内存

3. 选择合适的算法

python

from scipy import linalg
import numpy as np

# 对于对称正定矩阵，使用 Cholesky 分解更高效
n = 1000
A = np.random.random((n, n))
A = A @ A.T  # 创建对称正定矩阵
b = np.random.random(n)

# 通用 LU 分解
start_time = time.time()
x1 = linalg.solve(A, b)
lu_time = time.time() - start_time

# Cholesky 分解（对称正定矩阵）
start_time = time.time()
x2 = linalg.solve(A, b, assume_a='pos')
chol_time = time.time() - start_time

print(f"LU 分解耗时: {lu_time:.4f} 秒")
print(f"Cholesky 分解耗时: {chol_time:.4f} 秒")
print(f"性能提升: {lu_time / chol_time:.1f} 倍")
print(f"解的差异: {np.max(np.abs(x1 - x2)):.2e}")

实际应用示例

示例 1：数据预处理

python

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 生成带噪声的数据
np.random.seed(42)
n_samples = 1000
true_signal = np.sin(np.linspace(0, 4*np.pi, n_samples))
noise = np.random.normal(0, 0.3, n_samples)
noisy_data = true_signal + noise

# 数据清洗：移除异常值
z_scores = np.abs(stats.zscore(noisy_data))
threshold = 3
clean_data = noisy_data[z_scores < threshold]
clean_indices = np.where(z_scores < threshold)[0]

print(f"原始数据点数: {len(noisy_data)}")
print(f"清洗后数据点数: {len(clean_data)}")
print(f"移除的异常值: {len(noisy_data) - len(clean_data)}")

# 可视化结果
plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.plot(true_signal, label='真实信号', linewidth=2)
plt.title('真实信号')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 2)
plt.plot(noisy_data, alpha=0.7, label='带噪声数据')
plt.title('带噪声数据')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 3)
plt.plot(clean_indices, clean_data, 'g.', alpha=0.7, label='清洗后数据')
plt.title('清洗后数据')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 4)
plt.hist(noisy_data, bins=50, alpha=0.7, label='原始数据', density=True)
plt.hist(clean_data, bins=50, alpha=0.7, label='清洗后数据', density=True)
plt.title('数据分布对比')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

示例 2：简单的机器学习流程

python

import numpy as np
from scipy import stats
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

# 生成合成数据
np.random.seed(42)
n_samples = 200
X = np.random.uniform(-3, 3, n_samples)
y = 2 * X + 1 + np.random.normal(0, 0.5, n_samples)  # y = 2x + 1 + noise

# 数据预处理
X = X.reshape(-1, 1)  # sklearn 需要 2D 数组

# 分割数据
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 训练模型
model = LinearRegression()
model.fit(X_train, y_train)

# 预测
y_pred = model.predict(X_test)

# 评估模型
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)

print(f"模型参数:")
print(f"  斜率: {model.coef_[0]:.3f} (真实值: 2.000)")
print(f"  截距: {model.intercept_:.3f} (真实值: 1.000)")
print(f"模型性能:")
print(f"  均方误差: {mse:.3f}")
print(f"  R² 分数: {r2:.3f}")

# 统计检验
residuals = y_test - y_pred
_, p_value = stats.normaltest(residuals)
print(f"残差正态性检验 p 值: {p_value:.3f}")

# 可视化
plt.figure(figsize=(12, 4))

plt.subplot(1, 3, 1)
plt.scatter(X_train, y_train, alpha=0.6, label='训练数据')
plt.scatter(X_test, y_test, alpha=0.6, label='测试数据')
X_line = np.linspace(X.min(), X.max(), 100).reshape(-1, 1)
y_line = model.predict(X_line)
plt.plot(X_line, y_line, 'r-', linewidth=2, label='拟合直线')
plt.xlabel('X')
plt.ylabel('y')
plt.title('线性回归拟合')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(1, 3, 2)
plt.scatter(y_test, y_pred, alpha=0.6)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', linewidth=2)
plt.xlabel('真实值')
plt.ylabel('预测值')
plt.title('预测 vs 真实值')
plt.grid(True, alpha=0.3)

plt.subplot(1, 3, 3)
plt.hist(residuals, bins=20, density=True, alpha=0.7)
x_norm = np.linspace(residuals.min(), residuals.max(), 100)
y_norm = stats.norm.pdf(x_norm, np.mean(residuals), np.std(residuals))
plt.plot(x_norm, y_norm, 'r-', linewidth=2, label='正态分布')
plt.xlabel('残差')
plt.ylabel('密度')
plt.title('残差分布')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

最佳实践

1. 代码组织

python

# 推荐的导入方式
import numpy as np
import scipy as sp
from scipy import stats, optimize, integrate, linalg
import matplotlib.pyplot as plt

# 避免使用 from scipy import *
# 这会污染命名空间

# 使用有意义的变量名
data_samples = np.random.normal(0, 1, 1000)  # 好
x = np.random.normal(0, 1, 1000)  # 不够清晰

# 添加适当的注释
def calculate_confidence_interval(data, confidence=0.95):
    """
    计算数据的置信区间
    
    参数:
        data: 数据数组
        confidence: 置信水平 (默认 0.95)
    
    返回:
        (下界, 上界): 置信区间的下界和上界
    """
    alpha = 1 - confidence
    n = len(data)
    mean = np.mean(data)
    std_err = stats.sem(data)  # 标准误差
    t_value = stats.t.ppf(1 - alpha/2, n-1)
    margin_error = t_value * std_err
    
    return mean - margin_error, mean + margin_error

2. 错误处理

python

def safe_divide(a, b):
    """
    安全的除法操作，处理除零错误
    """
    try:
        result = np.divide(a, b)
        # 检查是否有无穷大或 NaN
        if np.any(np.isinf(result)) or np.any(np.isnan(result)):
            print("警告: 结果包含无穷大或 NaN 值")
        return result
    except ZeroDivisionError:
        print("错误: 除零操作")
        return None
    except Exception as e:
        print(f"未知错误: {e}")
        return None

# 使用示例
a = np.array([1, 2, 3, 4])
b = np.array([2, 0, 1, 2])
result = safe_divide(a, b)
print(f"安全除法结果: {result}")

3. 性能监控

python

import time
from functools import wraps

def timing_decorator(func):
    """
    装饰器：测量函数执行时间
    """
    @wraps(func)
    def wrapper(*args, **kwargs):
        start_time = time.time()
        result = func(*args, **kwargs)
        end_time = time.time()
        print(f"{func.__name__} 执行时间: {end_time - start_time:.4f} 秒")
        return result
    return wrapper

@timing_decorator
def matrix_multiplication(n):
    """
    矩阵乘法性能测试
    """
    A = np.random.random((n, n))
    B = np.random.random((n, n))
    return A @ B

# 测试不同大小的矩阵
for size in [100, 500, 1000]:
    print(f"\n矩阵大小: {size}x{size}")
    result = matrix_multiplication(size)

小结

在本章中，我们学习了 SciPy 的基础概念：

核心理念：基于 NumPy 数组、模块化设计、函数式编程
数据类型：标量、向量、矩阵、张量
数组操作：索引、切片、形状操作、拼接分割
函数模式：统计、优化、积分等常用模式
错误处理：常见错误类型和调试技巧
性能优化：向量化、内存管理、算法选择
最佳实践：代码组织、错误处理、性能监控

掌握这些基础概念后，您就可以开始学习 SciPy 的具体功能模块了。接下来，我们将在数组和矩阵操作中学习更高级的数组处理技术。

练习题

数组操作练习：创建一个 5x5 的随机矩阵，计算其行和、列和、对角线和
性能比较：比较使用循环和向量化操作计算数组元素平方的性能差异
错误处理：编写一个函数，安全地计算数组的对数，处理负数和零值
统计分析：生成两组随机数据，进行 t 检验比较它们的均值差异

python

# 练习题参考答案

# 1. 数组操作练习
import numpy as np

matrix = np.random.random((5, 5))
print(f"矩阵:\n{matrix}")
print(f"行和: {np.sum(matrix, axis=1)}")
print(f"列和: {np.sum(matrix, axis=0)}")
print(f"主对角线和: {np.trace(matrix)}")
print(f"反对角线和: {np.sum(np.fliplr(matrix).diagonal())}")

# 2. 性能比较
import time

n = 100000
data = np.random.random(n)

# 循环方式
start = time.time()
squared_loop = [x**2 for x in data]
loop_time = time.time() - start

# 向量化方式
start = time.time()
squared_vectorized = data**2
vectorized_time = time.time() - start

print(f"循环耗时: {loop_time:.4f}s")
print(f"向量化耗时: {vectorized_time:.4f}s")
print(f"性能提升: {loop_time/vectorized_time:.1f}倍")

# 3. 错误处理
def safe_log(arr):
    """
    安全计算对数，处理负数和零值
    """
    result = np.full_like(arr, np.nan, dtype=float)
    positive_mask = arr > 0
    result[positive_mask] = np.log(arr[positive_mask])
    
    if np.any(arr <= 0):
        print(f"警告: 发现 {np.sum(arr <= 0)} 个非正数值")
    
    return result

# 测试
test_data = np.array([-1, 0, 1, 2.718, 10])
log_result = safe_log(test_data)
print(f"安全对数结果: {log_result}")

# 4. 统计分析
from scipy import stats

np.random.seed(42)
group1 = np.random.normal(100, 15, 50)
group2 = np.random.normal(105, 15, 50)

t_stat, p_value = stats.ttest_ind(group1, group2)
print(f"两样本 t 检验:")
print(f"t 统计量: {t_stat:.3f}")
print(f"p 值: {p_value:.3f}")
print(f"结论: {'存在显著差异' if p_value < 0.05 else '无显著差异'}")

SciPy 基础概念 ​

SciPy 的核心理念 ​

1. 基于 NumPy 数组 ​

2. 模块化设计 ​

3. 函数式编程风格 ​

常用数据类型 ​

1. 标量（Scalars） ​

2. 一维数组（Vectors） ​

3. 二维数组（Matrices） ​

4. 多维数组（Tensors） ​

基本数组操作 ​

1. 数组索引和切片 ​

2. 数组形状操作 ​

3. 数组拼接和分割 ​

SciPy 常用函数模式 ​

1. 统计函数模式 ​

2. 优化函数模式 ​

3. 积分函数模式 ​

错误处理和调试 ​

1. 常见错误类型 ​

2. 调试技巧 ​

性能优化基础 ​

1. 向量化操作 ​

2. 内存管理 ​

3. 选择合适的算法 ​

实际应用示例 ​

示例 1：数据预处理 ​

示例 2：简单的机器学习流程 ​

最佳实践 ​

1. 代码组织 ​

2. 错误处理 ​

3. 性能监控 ​

小结 ​

练习题 ​

SciPy 基础概念

SciPy 的核心理念

1. 基于 NumPy 数组

2. 模块化设计

3. 函数式编程风格

常用数据类型

1. 标量（Scalars）

2. 一维数组（Vectors）

3. 二维数组（Matrices）

4. 多维数组（Tensors）

基本数组操作

1. 数组索引和切片

2. 数组形状操作

3. 数组拼接和分割

SciPy 常用函数模式

1. 统计函数模式

2. 优化函数模式

3. 积分函数模式

错误处理和调试

1. 常见错误类型

2. 调试技巧

性能优化基础

1. 向量化操作

2. 内存管理

3. 选择合适的算法

实际应用示例

示例 1：数据预处理

示例 2：简单的机器学习流程

最佳实践

1. 代码组织

2. 错误处理

3. 性能监控

小结

练习题