如何将向量分割为n个相似段？求R/Python可用的相关算法

给定一个含m个实数的向量，需将其分割为位置相邻且数值相似的段。相似性可定义为最小化段内数值变异性（如方差），例如向量[4, 4.2, 4, 18, 1, 2, 0.98, 15, 17]分割为4段时，结果为{[4, 4.2, 4], [18], [1, 2, 0.98],[15, 17]}。以下针对两个问题给出解决方案：

问题1：指定段数n的最优分割算法 #

**动态规划（Dynamic Programming）**是解决这类问题的标准方法，可精准找到最小化各段方差之和的最优分割方式。

算法逻辑 #

1. 预处理代价矩阵：提前计算所有可能子段的方差，存储在二维数组cost[i][j]中，代表从第i个元素到第j个元素组成的子段的方差（变异性代价）。
2. 动态规划状态定义：设dp[k][j]表示前j个元素分割为k段时的最小总代价。
3. 状态转移：dp[k][j] = min(dp[k-1][i] + cost[i+1][j])，其中i的范围为k-1 ≤ i < j（保证前k-1段至少有k-1个元素）。
4. 回溯分割点：计算完dp数组后，从dp[n][m]反向推导得到所有分割位置。

Python实现示例 #

import numpy as np

def segment_fixed_n(arr, n):
    m = len(arr)
    # 预处理所有子段的方差
    cost = np.zeros((m, m))
    for i in range(m):
        for j in range(i, m):
            sub = arr[i:j+1]
            cost[i][j] = np.var(sub)
    
    # 初始化DP数组
    dp = np.full((n+1, m), np.inf)
    # 分割为1段时的代价
    for j in range(m):
        dp[1][j] = cost[0][j]
    
    # 填充DP数组
    for k in range(2, n+1):
        for j in range(k-1, m):
            for i in range(k-2, j):
                if dp[k-1][i] + cost[i+1][j] < dp[k][j]:
                    dp[k][j] = dp[k-1][i] + cost[i+1][j]
    
    # 回溯找分割点
    splits = []
    current = m-1
    for k in range(n, 1, -1):
        for i in range(current-1, k-3, -1):
            if dp[k][current] == dp[k-1][i] + cost[i+1][current]:
                splits.append(i+1)
                current = i
                break
    splits.reverse()
    
    # 生成分段结果
    segments = []
    prev = 0
    for s in splits:
        segments.append(arr[prev:s])
        prev = s
    segments.append(arr[prev:])
    return segments

# 测试示例
arr = [4, 4.2, 4, 18, 1, 2, 0.98, 15, 17]
print(segment_fixed_n(arr, 4))

R实现示例 #

segment_fixed_n <- function(arr, n) {
  m <- length(arr)
  # 预处理子段方差
  cost <- matrix(0, nrow = m, ncol = m)
  for (i in 1:m) {
    for (j in i:m) {
      sub <- arr[i:j]
      cost[i,j] <- var(sub)
    }
  }
  
  # 初始化DP矩阵
  dp <- matrix(Inf, nrow = n+1, ncol = m)
  # 1段的情况
  dp[1, ] <- cost[1, ]
  
  # 填充DP
  for (k in 2:(n+1)) {
    for (j in k:m) {
      for (i in (k-1):(j-1)) {
        if (dp[k-1, i] + cost[i+1, j] < dp[k, j]) {
          dp[k, j] <- dp[k-1, i] + cost[i+1, j]
        }
      }
    }
  }
  
  # 回溯分割点
  splits <- c()
  current <- m
  for (k in (n+1):2) {
    for (i in (current-1):(k-1)) {
      if (dp[k, current] == dp[k-1, i] + cost[i+1, current]) {
        splits <- c(i, splits)
        current <- i
        break
      }
    }
  }
  
  # 生成分段
  segments <- list()
  prev <- 1
  for (s in splits) {
    segments[[length(segments)+1]] <- arr[prev:s]
    prev <- s+1
  }
  segments[[length(segments)+1]] <- arr[prev:m]
  return(segments)
}

# 测试
arr <- c(4, 4.2, 4, 18, 1, 2, 0.98, 15, 17)
print(segment_fixed_n(arr, 4))

问题2：自动确定最优段数的算法 #

这类问题需引入正则化惩罚项，避免段数过多导致的过度拟合（如m段时总变异性为0，但无实际意义）。以下是两种常用方案：

1. 带BIC正则化的动态规划 #

算法逻辑

目标函数结合总变异性与段数惩罚：总代价 = m*log(总方差和/m) + n*log(m)（BIC形式），遍历1到m的所有可能段数n，用动态规划计算对应最小总代价，选择代价最小的n作为最优段数。

2. PELT变化点检测算法 #

**Pruned Exact Linear Time（PELT）**是高效的序列变化点检测算法，通过剪枝策略优化动态规划，可自动找到最优分割点与段数，适合长序列处理。

Python实现示例（基于ruptures库）

import ruptures as rpt

def segment_auto(arr):
    # 采用L2模型（最小化方差和），PELT算法检测变化点
    model = "l2"
    algo = rpt.Pelt(model=model).fit(arr)
    # pen为惩罚系数，需根据数据调整
    breaks = algo.predict(pen=10)
    # 生成分段结果
    segments = []
    prev = 0
    for b in breaks:
        segments.append(arr[prev:b])
        prev = b
    return segments

# 测试示例
arr = [4, 4.2, 4, 18, 1, 2, 0.98, 15, 17]
print(segment_auto(arr))

注：需先安装第三方库pip install ruptures

R实现示例（基于changepoint库）

library(changepoint)

segment_auto <- function(arr) {
  # 使用PELT算法检测方差变化点，采用BIC惩罚
  cpt_result <- cpt.var(arr, method="PELT", penalty="BIC")
  breaks <- cpts(cpt_result)
  # 生成分段结果
  segments <- list()
  prev <- 1
  for (b in breaks) {
    segments[[length(segments)+1]] <- arr[prev:b]
    prev <- b+1
  }
  segments[[length(segments)+1]] <- arr[prev:length(arr)]
  return(segments)
}

# 测试
arr <- c(4, 4.2, 4, 18, 1, 2, 0.98, 15, 17)
print(segment_auto(arr))

注：需先安装第三方库install.packages("changepoint")

内容的提问来源于stack exchange，提问作者oska boska