布丰投针问题中长针（l>d）恰好k次相交的概率推导校验、模拟偏差排查及可视化咨询

嗨，我帮你梳理下推导、模拟里的问题，再聊聊可视化的思路：

一、推导逻辑的校验 #

先确认你的核心推导框架是没问题的：因为用了水平平行线，把原推导里的余弦换成正弦完全合理——只要把θ定义为针与水平线的夹角，这个转换不会改变推导的本质。

关于恰好k次相交的概率计算：
你设定$s \sim \mathcal{U}(0,1)$，令$r = l/d$（针长与间距的比值），恰好k次相交时θ的有效范围是：
$$\arcsin\left(\frac{k-s}{r}\right) \leq \theta < \arcsin\left(\frac{(k+1)-s}{r}\right)$$
对应的概率积分式（θ的概率密度为$2/\pi$，因为$\theta \sim \mathcal{U}(0,\pi/2)$）是：
$$P(k) = \int_0^1 \int_{\arcsin\left(\frac{k-s}{r}\right)}^{\arcsin\left(\frac{(k+1)-s}{r}\right)} \frac{2}{\pi} d\theta ds$$

但要注意定义域的有效性，这很可能是你积分求值出错的关键：

• 当$\frac{k-s}{r} > 1$或$\frac{k-s}{r} < -1$时，$\arcsin$无意义，所以积分的s范围其实不是全程0到1，而是要先筛选出s使得$\frac{k-s}{r} \in [-1,1]$且$\frac{k+1-s}{r} \in [-1,1]$
• 对于边界情况：
  • $k=0$时，θ的下限取0，积分范围改为$\int_0^1 \int_{0}^{\arcsin\left(\frac{1-s}{r}\right)} \frac{2}{\pi} d\theta ds$
  • $k=m$（$m=\lfloor r \rfloor$）时，θ的上限取$\pi/2$，积分范围改为$\int_0^1 \int_{\arcsin\left(\frac{m-s}{r}\right)}^{\pi/2} \frac{2}{\pi} d\theta ds$

二、模拟代码的偏差排查 #

你的模拟代码里有几个可能导致结果偏差的点，建议逐一排查：

1. m值定义错误：你写的m = math.ceil(r + 1)是错的，$m$应该是$\lfloor r \rfloor$（即math.floor(r)），比如$r=3$时$m=3$，错误的m会导致crossings数组长度冗余，统计结果失真。
2. 相交计数逻辑验证：可以加个手动测试用例：固定$y=0.5$，$\theta=\pi/2$（$\sin\theta=1$），$r=3$，此时针的x范围是0.5到3.5，会穿过x=1、2、3三条线，即3次相交，运行代码看是否返回no_of_crossings=3，如果不对，说明计数逻辑的条件判断有问题。
3. 理论值公式错误：你注释掉的theoretical_result(r)函数里的公式可能和积分结果不匹配，建议用符号积分工具（比如SymPy）先算出精确的理论表达式，再转换成Python代码，避免手动推导公式的错误。

以下是你的模拟代码（已修正m值的定义）：

import matplotlib.pyplot as plt
import numpy as np
import random as random
import math

no_of_needles = 10000
r = 3
m = math.floor(r)  # 修正m的定义
crossings = np.zeros(m + 1)

plt.figure(1)
# 绘制间距为1的平行竖线
for i in range(0, m+2):
    plt.axvline(i, color='black')

def generate_random_colors(m):
    # 生成m个随机RGB颜色
    colors = [(random.random(), random.random(), random.random()) for _ in range(m)]
    return colors

colors = generate_random_colors(len(crossings))
ys = np.arange(1, r+1, 1.0)

for i in range(no_of_needles):
    # 随机位置y~U(0,1)
    y = random.uniform(0, 1)
    # 随机角度θ~U(0, π/2)
    theta = random.uniform(0, np.pi/2)
    x_coords = [y, y + r * np.sin(theta)]
    y_coords = [0, r * np.cos(theta)]

    # 计算相交次数
    y_ns = np.insert(ys, 0, y)
    theta_ns = np.arcsin((y_ns - y) / r)
    theta_ns_new = theta_ns.searchsorted(theta)
    no_of_crossings = theta_ns_new - 1

    crossings[no_of_crossings] += 1
    plt.plot(x_coords, y_coords, color=colors[no_of_crossings], label=no_of_crossings)

def legend_without_duplicate_labels(figure):
    handles, labels = plt.gca().get_legend_handles_labels()
    by_label = dict(zip(labels, handles))
    figure.legend(by_label.values(), by_label.keys(), loc='lower right')

legend_without_duplicate_labels(plt.gcf())
plt.title('Needle orientation and positions colored by no of crossings')

plt.figure(2)
# 绘制相交次数的概率分布
plt.scatter(np.arange(0, len(crossings)), crossings/no_of_needles, color='black', marker='o')
plt.plot(np.arange(0, len(crossings)), crossings/no_of_needles, color='red')
plt.title(r'Histogram of crossings')
plt.grid(True)
plt.show()

三、可视化思路 #

要可视化恰好k次相交的积分约束，推荐绘制二维样本空间图：

• 横轴：s（0到1，对应针的位置偏移）
• 纵轴：θ（0到π/2，对应针的角度）
• 对于每个k，画出两条曲线：$\theta = \arcsin\left(\frac{k-s}{r}\right)$ 和 $\theta = \arcsin\left(\frac{(k+1)-s}{r}\right)$，两条曲线之间的区域（加上s的有效范围）就是恰好k次相交的样本空间

举个例子，当$r=3$、$k=1$时：

• 曲线1：$\theta = \arcsin\left(\frac{1-s}{3}\right)$，s从0到1时，θ从$\arcsin(1/3)≈0.3398$到0
• 曲线2：$\theta = \arcsin\left(\frac{2-s}{3}\right)$，s从0到1时，θ从$\arcsin(2/3)≈0.7297$到$\arcsin(1/3)≈0.3398$
• 用Matplotlib的fill_between函数填充这两条曲线之间的区域，就能直观看到恰好1次相交的样本空间占比

对于边界情况：

• $k=0$时，填充θ从0到$\arcsin\left(\frac{1-s}{3}\right)$的区域
• $k=3$时，填充θ从$\arcsin\left(\frac{3-s}{3}\right)$到$\pi/2$的区域

备注：内容来源于stack exchange，提问作者Tanamas