from scipy.stats import poisson
# On average, 1 block is found on average every
# "block period" which is 10 minutes
avg_blocks_per_period = 1
# but we actually found (for example) 5 blocks
obs_blocks = 5
prob = poisson.pmf(k=obs_blocks, mu=avg_blocks_per_period)
print(f"Probability of {obs_blocks} blocks in {avg_blocks_per_period} ten minute period is approximately {prob*100:.2f}%")
This should result in:
Probability of 5 blocks in 1 ten minute period is approximately 0.31%
At that rate, a 5-block run in a single 10 minute period should happen roughly a few times a difficulty period, even when the overall average is still just at 10 minutes.
I think this is correct anyway, it's been awhile since I had to do these kinds of stats.
The above doesn't work well if you want more fine grained stats, but you can, for example, see what the chances of 145 blocks a day are if you set 144 expected per day. There are other ways to get down to inter-block periods, but you have to adjust everything.
Probability of 5 blocks in 1 ten minute period is approximately 0.31%