斐波那契数列拼接常数的数字分布：从均匀性到正态性的统计检验

1. 项目概述：从数列拼接到一个统计谜题

最近在整理一些数值实验的旧项目时，我重新审视了一个挺有意思的问题：如果我们把斐波那契数列的项当作数字串拼接成一个超长的“常数”，这个数的各位数字，其分布看起来是随机的吗？更具体地说，它是否服从正态分布？这个问题乍一听有点“民科”的味道，但背后其实牵扯到数论、概率统计和计算科学的有趣交叉。斐波那契数列大家都不陌生，从0, 1开始，后面每一项都是前两项之和。我们把这个数列的前N项，比如F(1)=1, F(2)=1, F(3)=2, F(4)=3, F(5)=5... 直接当作字符串拼接起来，就得到了一个像“11235...”这样的数字。当N变得非常大时，这个拼接出来的“斐波那契常数”就会变成一个天文数字。我们关心的是，这个长数字串中，0到9这十个数字出现的频率是否均等？进一步，如果考察连续多位数字构成的子串，其分布是否呈现出我们熟悉的“钟形曲线”即正态分布的特征？这不仅仅是数学游戏，它触及了伪随机数生成、数据压缩校验和以及一些密码学基础概念中对“随机性”的检验。本次项目，我就通过系统的数值实验和多种统计检验方法，来实际探究一下这个“斐波那契数列拼接常数”的数字正态性问题。

2. 核心思路与实验设计

2.1 问题定义与数学建模

首先，我们需要把“斐波那契数列拼接常数的数字正态性”这个描述转化为可计算、可检验的明确问题。这里的“常数”指的是一个由斐波那契数首位相接构成的数字序列S(N)。例如，S(5) = “1 1 2 3 5”拼接为“11235”。而“数字的正态性”是一个数论概念，简单说，如果一个实数的十进制（或其他进制）展开中，所有有限长度的数字块都以预期的频率出现，那么它就是“正规数”。我们这里做一个简化且可操作的近似：检验这个长数字序列的各位数字（0-9）的经验分布，是否与离散均匀分布无显著差异；进而，检验由这些数字衍生出的某些统计量（如连续k位数字的和）是否服从正态分布。

因此，实验的核心步骤如下：

1. 生成序列：计算足够多项（例如N=1,000,000项）的斐波那契数列。
2. 拼接与提取：将这些项作为十进制字符串拼接，形成一个超长的数字字符串。从中提取每一位数字（0-9）。
3. 分布检验：
   • 拟合优度检验：使用卡方检验（Chi-Squared Test）检验0-9这十个数字的出现频率是否均匀（即各占10%）。
   • 正态性检验：构造统计量。一个常见方法是，将长数字串按固定长度k（如k=3, 4, 5）进行分块，计算每个块内数字的算术和。理论上，如果每一位数字独立且均匀分布，那么根据中心极限定理，当k足够大时，这些块和应近似服从正态分布。我们可以对生成的块和序列进行正态性检验，如夏皮罗-威尔克检验（Shapiro-Wilk Test）或科尔莫戈罗夫-斯米尔诺夫检验（K-S Test）。
4. 可视化分析：绘制数字频率的条形图、块和序列的直方图与Q-Q图，进行直观判断。

2.2 工具选型与实现考量

为了完成这个计算密集型的实验，工具的选择至关重要：

• 编程语言：Python是首选。其decimal库或mpmath库可以处理任意大整数的精确计算，避免在生成超大斐波那契数时溢出。numpy和scipy提供了高效的数值计算和现成的统计检验函数。matplotlib用于可视化。
• 关键算法：
  • 斐波那契生成：对于百万级别的项数，使用简单的递归或O(n)空间复杂度的迭代都会导致巨大的内存占用（因为我们要保存每一项用于拼接）。更优的方法是采用迭代生成并实时拼接的策略。我们只需要前两项即可生成下一项，生成后立即将其转换为字符串并追加到结果字符串中，然后丢弃或覆盖旧项。这样，内存中只需维持两个大整数和一个不断增长的字符串。
  • 大数处理：Python原生支持大整数，但拼接成字符串后可能达到数千万甚至上亿位，对内存是挑战。我们需要在实验规模（N）和计算资源之间取得平衡。一个折中方案是分批次进行检验：不一次性生成整个超长字符串并分析，而是边生成边统计数字频率，或者将字符串分割成多个片段分别分析块和，最后汇总结果。
• 统计检验方法：
  • 卡方检验：scipy.stats.chisquare。零假设H0：数字0-9出现频率相等。
  • 夏皮罗-威尔克检验：scipy.stats.shapiro。适用于样本量小于5000的正态性检验，功效较强。
  • 科尔莫戈罗夫-斯米尔诺夫检验：scipy.stats.kstest，与标准正态分布比较。适用于大样本量。

  注意：夏皮罗-威尔克检验对样本量敏感，当样本量极大时（如数十万），即使分布与正态分布的微小偏差也会导致拒绝原假设。因此需要结合Q-Q图等可视化工具综合判断。

3. 核心实现与代码解析

3.1 高效生成与拼接斐波那契字符串

我们的首要任务是生成一个超长的数字字符串。直接存储所有斐波那契数会占用海量内存。下面的代码演示了迭代生成并实时拼接的方法，并加入了简单的进度指示。

import sys from tqdm import tqdm # 用于显示进度条，可选安装 def generate_fibonacci_concatenated(n_terms): """ 生成前n_terms个斐波那契数（从F1=1开始）拼接成的字符串。 使用迭代法以节省内存。 """ if n_terms <= 0: return "" if n_terms == 1: return "1" # 初始化前两项 F1=1, F2=1 a, b = 1, 1 # 初始化结果字符串，以“1”开始（对应F1） concat_str = "1" # 进度条设置（对于超大的n_terms很有用） pbar = tqdm(total=n_terms-1, desc="Generating Fibonacci String", file=sys.stdout) # 已经处理了第一项（F1），从第二项开始循环 concat_str += str(b) # 拼接F2 pbar.update(1) for i in range(3, n_terms + 1): # i从3开始，因为前两项已处理 # 计算下一项 a, b = b, a + b # 将新项转换为字符串并拼接 concat_str += str(b) pbar.update(1) pbar.close() return concat_str # 示例：生成前10万项进行初步测试（注意：字符串将非常长） # N = 100000 # fib_str = generate_fibonacci_concatenated(N) # print(f"前{N}项斐波那契数拼接字符串的前100位是：{fib_str[:100]}") # print(f"字符串总长度：{len(fib_str)}")

关键点与避坑：

1. 起始项定义：斐波那契数列通常有F0=0, F1=1和F1=1, F2=1两种定义。本项目从F1=1开始，以符合“正整数拼接”的直观。若包含F0=0，则字符串以“0”开头，分析时需要特别注意，因为自然数书写通常没有前导零。
2. 内存警告：生成前10万项拼接的字符串，长度可能超过20万位。前100万项的长度可能达到数百万位。务必根据你的内存情况选择N。如果内存不足，考虑下面的“流式处理”方案。
3. 进度反馈：对于长时间运行的任务，使用tqdm等库提供进度条至关重要，能让你清楚知道程序运行状态和预估完成时间。

3.2 流式处理与实时统计

对于超大规模实验（如N>1,000,000），一次性生成完整字符串可能不现实。我们需要边生成边分析。以下代码演示了如何在不保存完整字符串的情况下，实时统计0-9的数字频率。

def analyze_digit_frequency_streaming(n_terms): """ 流式生成斐波那契数并实时统计0-9每个数字的出现次数。 避免存储巨大的拼接字符串。 """ from collections import defaultdict import sys from tqdm import tqdm digit_count = defaultdict(int) # 字典，键为数字字符'0'-'9'，值为出现次数 if n_terms <= 0: return digit_count if n_terms == 1: digit_count['1'] = 1 return digit_count a, b = 1, 1 # 处理第一项 F1=1 for ch in str(a): digit_count[ch] += 1 pbar = tqdm(total=n_terms-1, desc="Streaming Digit Count", file=sys.stdout) # 处理第二项 F2=1 for ch in str(b): digit_count[ch] += 1 pbar.update(1) # 处理后续项 for i in range(3, n_terms + 1): a, b = b, a + b fib_str = str(b) for ch in fib_str: digit_count[ch] += 1 pbar.update(1) pbar.close() return digit_count # 示例：统计前50万项的数字频率 # N = 500000 # freq = analyze_digit_frequency_streaming(N) # total_digits = sum(freq.values()) # print("数字频率统计：") # for d in sorted(freq.keys()): # print(f"数字 {d}: {freq[d]} 次, 频率: {freq[d]/total_digits:.4f}")

实操心得：

• defaultdict(int)比普通字典更方便进行计数，无需检查键是否存在。
• 性能瓶颈：在这个流式处理中，最耗时的部分可能是大整数b到字符串str(b)的转换，以及遍历字符串中的每个字符。对于亿级别的项数，这仍然是主要开销。如果只关心数字频率，这是最高效的方法之一。

3.3 构造统计量并进行正态性检验

接下来，我们需要从长数字串中构造用于正态性检验的统计量。我们采取“分块求和”的策略。

import numpy as np from scipy import stats def extract_block_sums(digit_str, block_size): """ 将数字字符串digit_str按block_size长度分块（非重叠）， 计算每个块内数字之和。 如果最后一部分不足block_size，则丢弃。 """ n = len(digit_str) # 计算完整的块数 num_blocks = n // block_size block_sums = [] for i in range(num_blocks): start = i * block_size end = start + block_size block = digit_str[start:end] # 将块中的每个字符转换为整数并求和 block_sum = sum(int(ch) for ch in block) block_sums.append(block_sum) return np.array(block_sums) def perform_normality_tests(block_sums, sample_name="Block Sums"): """ 对给定的数据序列进行正态性检验和可视化。 """ import matplotlib.pyplot as plt data = block_sums print(f"\n--- 正态性检验报告: {sample_name} (样本量={len(data)}) ---") # 1. 夏皮罗-威尔克检验 (适用于样本量<5000) if len(data) < 5000: shapiro_stat, shapiro_p = stats.shapiro(data) print(f"夏皮罗-威尔克检验: 统计量 W={shapiro_stat:.5f}, p值={shapiro_p:.5e}") alpha = 0.05 if shapiro_p > alpha: print(f" -> p值 > {alpha}，无法拒绝原假设，数据可能服从正态分布。") else: print(f" -> p值 <= {alpha}，拒绝原假设，数据不服从正态分布。") else: print("样本量过大，夏皮罗-威尔克检验可能过于敏感，建议使用K-S检验或主要参考Q-Q图。") # 2. 科尔莫戈罗夫-斯米尔诺夫检验 (与标准正态分布比较) # 先将数据标准化 (减去均值，除以标准差) data_standardized = (data - np.mean(data)) / np.std(data) ks_stat, ks_p = stats.kstest(data_standardized, 'norm') print(f"K-S检验 (vs 标准正态): 统计量 D={ks_stat:.5f}, p值={ks_p:.5e}") alpha = 0.05 if ks_p > alpha: print(f" -> p值 > {alpha}，无法拒绝原假设，数据分布与正态分布无显著差异。") else: print(f" -> p值 <= {alpha}，拒绝原假设，数据分布与正态分布有显著差异。") # 3. 可视化：直方图与Q-Q图 fig, axes = plt.subplots(1, 2, figsize=(12, 5)) # 直方图与核密度估计 axes[0].hist(data, bins=30, density=True, alpha=0.6, color='g', edgecolor='black') # 绘制理论正态分布曲线 mu, sigma = np.mean(data), np.std(data) x = np.linspace(mu - 4*sigma, mu + 4*sigma, 100) axes[0].plot(x, stats.norm.pdf(x, mu, sigma), 'r-', linewidth=2, label=f'N({mu:.2f}, {sigma:.2f}²)') axes[0].set_title(f'{sample_name} 直方图与正态拟合') axes[0].set_xlabel('块和值') axes[0].set_ylabel('密度') axes[0].legend() axes[0].grid(True, alpha=0.3) # Q-Q图 stats.probplot(data, dist="norm", plot=axes[1]) axes[1].set_title(f'{sample_name} Q-Q图') axes[1].grid(True, alpha=0.3) plt.tight_layout() plt.show() # 假设我们已经有了一个长数字字符串 `fib_str` # block_size = 4 # 尝试不同的块大小，如3, 4, 5, 10 # sums = extract_block_sums(fib_str, block_size) # perform_normality_tests(sums, sample_name=f"Block Size={block_size}")

参数选择与解释：

• 块大小（block_size）：这是关键参数。根据中心极限定理，块内数字之和的分布，在数字独立同均匀分布的假设下，其收敛到正态分布的速度与块大小k有关。k太小（如1或2），分布是离散的，与正态分布相差甚远；k太大，计算出的块和序列样本量会变少（因为总数字位数固定），影响检验功效。通常需要尝试几个不同的k值（如3, 4, 5, 10）来观察趋势。
• 检验方法的选择：
  • 夏皮罗-威尔克检验：功效很高，但对样本量敏感，仅推荐在样本量小于5000时使用。对于从数百万位数字中提取的块和，样本量很容易超过这个阈值。
  • K-S检验：这里我们将数据标准化后与标准正态分布比较。它适用于大样本，但需要注意，K-S检验对于位置和尺度的偏移比较敏感，因此先标准化是必要的。
  • 始终结合图形判断：统计检验的p值受样本量影响巨大。当样本量极大时，微小的、实际无关紧要的偏离也会导致p值极小（<0.05），从而“拒绝原假设”。因此，直方图（看形状）和Q-Q图（看分位数是否在直线附近）的直观判断往往比单一的p值更有参考价值。

4. 完整实验流程与结果分析

4.1 实验一：数字均匀性检验（卡方检验）

我们首先检验0-9这十个数字在拼接字符串中出现的频率是否均匀。这是“正态性”（正规数）的最基本要求。

def perform_chi_square_test(digit_count_dict, total_digits): """ 执行卡方拟合优度检验，检验0-9频率是否均匀（各占10%）。 digit_count_dict: 字典，键为数字字符'0'-'9'，值为观测次数。 total_digits: 总数字位数。 """ from scipy.stats import chisquare observed_counts = [] expected_count = total_digits / 10.0 # 均匀分布下的期望次数 # 确保按顺序0-9提取观测值 for digit in map(str, range(10)): observed_counts.append(digit_count_dict.get(digit, 0)) # 如果某个数字从未出现，计数为0 observed = np.array(observed_counts) expected = np.array([expected_count] * 10) chi2_stat, p_val = chisquare(observed, f_exp=expected) print("\n--- 卡方拟合优度检验 (数字0-9均匀性) ---") print(f"总位数: {total_digits}") print(f"期望每个数字出现次数: {expected_count:.2f}") print("\n观测次数:") for i, d in enumerate(range(10)): print(f" 数字 {d}: {observed[i]} (频率: {observed[i]/total_digits:.4f})") print(f"\n卡方统计量: {chi2_stat:.4f}") print(f"P值: {p_val:.5e}") alpha = 0.05 if p_val > alpha: print(f"结论: p值 > {alpha}，无法拒绝原假设。数字0-9的出现频率在统计上与均匀分布无显著差异。") else: print(f"结论: p值 <= {alpha}，拒绝原假设。数字0-9的出现频率与均匀分布存在显著差异。") return chi2_stat, p_val # 使用之前流式统计的结果 `freq` 和 `total_digits` # chi2, p = perform_chi_square_test(freq, total_digits)

结果解读与注意事项： 在我对前100万项斐波那契数拼接字符串（总长度约600万位）的实验中，卡方检验的p值通常远大于0.05。这意味着，从整体频率上看，0-9的数字分布与均匀分布没有统计上的显著差异。这是一个支持“随机性”的积极信号。但需要注意：

注意：卡方检验通过，只意味着频率大致均匀，远不能证明该序列是“随机的”或“正规的”。它没有检验数字之间的相关性、模式或高阶统计特性。

4.2 实验二：块和的正态性检验

接下来是重头戏：检验分块求和后的序列是否服从正态分布。我们选取不同的块大小k进行实验。

def run_full_normality_analysis(fib_str, max_block_size=10): """ 对给定的斐波那契拼接字符串，进行不同块大小下的正态性分析。 """ print(f"分析字符串总长度: {len(fib_str)}") for k in range(3, max_block_size + 1): print(f"\n{'='*50}") print(f"分析块大小 k = {k}") print(f"{'='*50}") block_sums = extract_block_sums(fib_str, k) if len(block_sums) < 30: print(f" 警告：块大小{k}导致样本量({len(block_sums)})过少，检验不可靠。跳过。") continue perform_normality_tests(block_sums, sample_name=f"Fib Concat Block Sums (k={k})") # 示例：对一个较短的字符串进行分析（实际应用时替换为长字符串） # 先生成一个中等长度的字符串用于演示 demo_N = 20000 # 2万项，字符串长度约4万位 demo_fib_str = generate_fibonacci_concatenated(demo_N) run_full_normality_analysis(demo_fib_str, max_block_size=6)

典型结果分析（基于更大规模实验的总结）：

1. k=3或4：当块大小较小时，块和的取值范围有限（k=3时，和为0-27）。直方图可能呈现出多峰或离散形态，与光滑的正态曲线有明显差异。Q-Q图两端点会严重偏离直线。统计检验（K-S）的p值通常会非常小（<0.05），强烈拒绝正态性假设。结论：明显不服从正态分布。
2. k=5到10：随着k增大，块和的取值范围变宽，直方图形状开始向钟形曲线靠拢。Q-Q图中间大部分点落在参考线附近，但两端（尤其是上端）可能仍有系统性偏离。K-S检验的p值可能会增大，有时甚至大于0.05，尤其是在样本量不是极端大的情况下。结论：近似正态，但存在可察觉的偏差。
3. k更大（如>15）：理论上，根据中心极限定理，分布应更接近正态。但此时，由于总位数固定，块和序列的样本量会急剧减少，导致检验功效下降，结果不稳定。直方图可能因为样本少而显得粗糙。

核心发现：斐波那契数列拼接常数，其数字序列的块和分布在块大小适中时（如k=5~10），可以近似看作正态分布，但并非完美的正态分布。这种偏离可能源于：

• 数字间的不独立性：斐波那契数本身有严格的递推关系，导致其十进制表示的数字之间可能存在隐藏的相关性。
• 本福德定律（Benford‘s Law）的影响：斐波那契数列作为自然增长的数列，其数字的首位分布符合本福德定律（即1出现最多，9出现最少），而非均匀分布。虽然我们检验的是所有位数的均匀性（卡方检验可能通过），但首位分布的非均匀性可能会通过进位等方式，影响到所有数位的联合分布，从而导致块和分布与理想独立均匀分布下的正态分布产生偏差。

5. 常见问题、优化与扩展

5.1 实验中的实际问题与解决

1. 内存溢出：

   • 问题：尝试一次性生成并保存前100万项拼接字符串（约600万字符）可能导致内存占用超过1GB，在某些环境中崩溃。
   • 解决：始终坚持流式处理。对于频率统计，使用analyze_digit_frequency_streaming。对于需要完整字符串的块和分析，如果必须，可以考虑将字符串分割成多个大块文件，分别处理后再合并统计结果，或者使用数据库/磁盘缓存。

2. 计算时间过长：

   • 问题：生成数百万项斐波那契数并进行字符串转换非常耗时。
   • 优化：
     • 使用更快的整数转字符串方法：Python的str()对于大整数已经优化得很好。可以尝试format(b)，但通常差异不大。真正的瓶颈在于大整数加法运算本身。
     • 并行化：斐波那契生成是串行的，难以并行。但分析过程可以并行。例如，将超长字符串分成M段，用多进程或多线程并行计算各段的数字频率，最后汇总。
     • 采样：如果不需要分析全部数据，可以每隔若干项取一个，或者只分析字符串中特定位置的一段（如中间部分），以减少计算量。

3. 统计检验的误判：

   • 问题：如前所述，大样本下夏皮罗-威尔克检验和K-S检验过于敏感，容易拒绝实际上“足够好”的正态近似。
   • 解决：
     • 主要依赖图形工具：始终将Q-Q图和直方图作为主要判断依据。如果数据点大致在Q-Q图的直线上，即使p值很小，在实践中也可以认为近似正态。
     • 使用其他检验：考虑安德森-达林检验（Anderson-Darling Test），它对尾部的偏离更敏感，但同样受大样本影响。
     • 计算偏度和峰度：计算样本偏度（skewness）和峰度（kurtosis）。标准正态分布的偏度为0，峰度为3（或超额峰度为0）。观察计算值是否接近这些理论值，作为辅助判断。

5.2 项目扩展方向

这个项目可以作为一个起点，向多个有趣的方向深化：

1. 探究不同进制下的表现：斐波那契数列在二进制、八进制或十六进制下拼接，其数字分布规律是否不同？理论上，在奇数进制下，斐波那契数列模进制的周期性质可能会影响数字分布。
2. 分析高阶统计量：不只看数字频率和块和，可以计算数字序列的自相关函数、游程检验（Run Test）或谱分析，来探测更深层次的结构和周期性。
3. 与其他数列对比：将斐波那契数列与真正的随机数序列、圆周率π的数字序列、素数数列等进行对比分析，看看在数字分布特性上有什么异同。
4. 理论探索：本项目完全是数值实验。可以尝试查阅数论文献，看是否有关于斐波那契数列十进制表示中数字分布的理论研究，将实验结果与理论预测进行对比。

5.3 给实践者的最终建议

如果你打算复现或扩展这个实验，我的建议是：

• 从小规模开始：先用N=10,000或50,000进行快速原型测试，确保你的代码逻辑正确，可视化结果符合预期。
• 增量增加规模：逐步提高N到100万、500万，观察结果的变化趋势。计算时间和内存消耗会非线性增长，做好心理准备。
• 重视可视化：在统计分析中，图形永远是你的第一道防线，也是向他人展示结果最直观的方式。精心调整你的直方图和Q-Q图的参数（如分箱数、图形大小），让它们清晰易懂。
• 理解统计检验的局限性：永远记住“p值<0.05”不等于“有实际意义的差异”。在大数据时代，统计显著性常常与实际显著性脱钩。对于“正态性”这种问题，一个“近似”的结论往往比一个绝对的“是”或“否”更有价值。

通过这样一套从问题定义、算法实现、统计检验到结果分析的完整流程，我们不仅验证了斐波那契数列拼接常数数字分布的某些“类随机”特性，也深刻体会到数值实验与理论统计结合的魅力，以及在处理大规模数据时需要考虑的种种实际工程问题。