我們知道，每次拔號(hào)上網(wǎng)運(yùn)營(yíng)商會(huì)給你分配一個(gè)IP地址，通常每次都不一樣。
我記錄了355 次拔號(hào)得到的IP地址。其中：

• 237 個(gè)地址是各出現(xiàn)了（僅）1次的；
• 50 個(gè)地址是各出現(xiàn)了（僅）2次的；
• 6 個(gè)地址是各出現(xiàn)了（僅）3次的。

因?yàn)楣W(wǎng)IP是稀缺資源，于是就想琢磨一下我這里運(yùn)營(yíng)商的IP地址池到底有多大。根據(jù)這些數(shù)據(jù)能否估算出地址池的大小N呢？

為了漂亮的封面

簡(jiǎn)單的數(shù)學(xué)模型

上面的實(shí)驗(yàn)中總共得到了293 (237 + 50 + 6)個(gè)IP地址，反過(guò)來(lái)，未被選中的地址是 N - 293個(gè)。
于是可知：

• 經(jīng)過(guò)355次拔號(hào)，某一特定地址不被選中的概率是 (N - 293)/N；
• 另一方面，僅拔號(hào)一次，某一特定的地址不被選中的概率是 (N - 1) / N。

于是我們可以得到下面的方程，并在Mathematica里求數(shù)值解（注意要限定實(shí)數(shù)域，否則會(huì)得到很多復(fù)數(shù)解）。

圖：簡(jiǎn)單模型的方程及其解

顯然 N > 293，于是有意義的解是 891.992。這是一個(gè)類似于平均值（數(shù)學(xué)期望）的數(shù)。
我們可以說(shuō)，N在892附近。至于有多近，現(xiàn)在還不清楚，N可能是800，也可能是1000。

不過(guò)，這個(gè)模型里沒(méi)有用到[237, 50, 6]這個(gè)1，2，3次的分布。這樣做對(duì)嗎？

簡(jiǎn)單的模擬實(shí)驗(yàn)

我準(zhǔn)備用Python模擬一下，用到下面一些package。

import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from multiprocessing import Pool

numpy里有一些現(xiàn)成的函數(shù)，用很短的代碼就可以完成模擬；pandas用來(lái)畫圖，最后為了提升速度會(huì)用到多進(jìn)程。代碼在Jupyter Notebook里運(yùn)行。

圖：簡(jiǎn)單的模擬實(shí)驗(yàn)

一次實(shí)驗(yàn)的代碼其實(shí)只有一行，流程是：

• 生成355個(gè)0~n之間的隨機(jī)數(shù)；
• 統(tǒng)計(jì)每個(gè)數(shù)出現(xiàn)的頻次；
• 統(tǒng)計(jì)每個(gè)頻次出現(xiàn)的頻次（即出現(xiàn)了0, 1, 2, 3... 次的數(shù)有多少個(gè)），其中出現(xiàn)0次的數(shù)目是不需要的。

把實(shí)驗(yàn)得到的頻次和目標(biāo)頻次進(jìn)行比較，相等就把計(jì)數(shù)加1. 對(duì)于每個(gè)n值，試驗(yàn)10000次，看看有多少命中的。

可以看到峰值確實(shí)是在900附近，大約是0.6%的命中率。也就是說(shuō)，對(duì)于N在900附近時(shí)，從中取355次地址，大約有0.6%的可能會(huì)出現(xiàn)[237, 50, 6]這樣的分布。

盡管這個(gè)值很低，但我們要看相對(duì)的高低，在900附近，是最有可能出現(xiàn)這個(gè)分布的。

程序的改進(jìn)

前面的代碼里可能會(huì)匹配到[237, 50, 6, 2]（即還有2個(gè)地址各出現(xiàn)了4次）這樣的分布，所以嚴(yán)謹(jǐn)一點(diǎn)，我們可以把匹配目標(biāo)設(shè)為[237, 50, 6, 0, 0, 0]。再長(zhǎng)也沒(méi)有必要了，因?yàn)楦哳l次的數(shù)量降低是很快的。

把程序改為多進(jìn)程模式的。4核8線程的CPU把Pool大小設(shè)為8,實(shí)測(cè)CPU占用率能達(dá)到100%（而之前確實(shí)是1/8左右的占用率）。把每一個(gè)N值的模擬次數(shù)增加到2千萬(wàn)，這樣可以得到較為平滑的曲線。

圖：改進(jìn)的模擬程序

這時(shí)可以得到最大概率對(duì)應(yīng)的N值是891, 和簡(jiǎn)單的數(shù)學(xué)模型非常一致。

概率公式和概率積分

根據(jù)“不盡相異元素的全排列”可以推出，出現(xiàn)[237, 50, 6]這個(gè)結(jié)果的概率p[n]，和n的關(guān)系，如下圖。

用這個(gè)結(jié)果畫出的曲線（圖略）和模擬的完全一致，但速度快得多，而且是完全平滑的。

直接把這個(gè)概率“積分”（嚴(yán)格講是累加），得到的結(jié)果是1.6. 為什么不是1呢？

因?yàn)檫@個(gè)概率p(n)表示的是對(duì)于某一個(gè)n，出現(xiàn)[237, 50, 6]這個(gè)結(jié)果的概率（它的反面是出現(xiàn)其它頻次的概率）；而我們需要的概率q(n)是指，[237, 50, 6]這個(gè)結(jié)果對(duì)應(yīng)的N，100%會(huì)落在比如說(shuō)293~10000之間，那么它落在n的概率是多大呢？

讓我們假定這兩個(gè)概率是正相關(guān)的，這樣把p(n)的積分歸一化，就可以得到q(n)。于是我們用概率公式得到下面歸一化的積分曲線。

圖：概率公式和概率積分

從圖中可以估算出，N 98%的可能會(huì)落在700~1200之間。曲線的陡峭程度反映了N的集中程度（如果取樣次數(shù)比355次更多，結(jié)果應(yīng)該會(huì)更集中）。

伯努利試驗(yàn)和二項(xiàng)分布

其實(shí)這個(gè)問(wèn)題是概率論中的一類典型問(wèn)題，可以和伯努利試驗(yàn)及二項(xiàng)分布對(duì)應(yīng)上。

在每次拔號(hào)中，選中某一地址為成功，否則為失敗。m次拔號(hào)中，某個(gè)地址被選中k次的分布符合二項(xiàng)分布 b(k, m, p)，其中p是被選中的概率（對(duì)于大小為n的地址池，p為1/n）。

于是，m次拔號(hào)中，被選中k次的地址的總數(shù)的期望值為：b(k, m, p) * n。具體如下圖所示：

圖：基于二項(xiàng)分布的計(jì)算

上圖中給出了出現(xiàn)1~5次的地址總數(shù)隨地址池大小n的分布。下圖則給出了目標(biāo)頻次[237, 50, 6]附近（+/-10%）所對(duì)應(yīng)的n值。比如，當(dāng)n取值700~1150時(shí)，最有可能出現(xiàn)237個(gè)1次地址。因?yàn)閳D中是期望值，所以并不意味著在這個(gè)范圍之外不會(huì)出現(xiàn)237個(gè)1次地址，只是圖中區(qū)域是概率最高的部分。這也可以理解，為什么根據(jù)3個(gè)條件得出的范圍并不是嚴(yán)格一致。

圖：基于二項(xiàng)分布的結(jié)果

結(jié)語(yǔ)

本文對(duì)“我這里”的公網(wǎng)IP地址的數(shù)量這一蛋疼的問(wèn)題，從數(shù)學(xué)上進(jìn)行了推演，從程序上進(jìn)行了模擬。數(shù)學(xué)簡(jiǎn)潔但高深，模擬笨拙但直觀。殊途同歸，兩者得出了基本一致的結(jié)果。

本題可以作為碼農(nóng)編程的一個(gè)小練習(xí)，比如有人用了PHP語(yǔ)言，因?yàn)樗懈呔鹊臄?shù)值計(jì)算；有人用了R語(yǔ)言，因?yàn)樗鞋F(xiàn)成的統(tǒng)計(jì)函數(shù)；繼續(xù)的改進(jìn)甚至可以考慮使用GPU提升模擬速度。

至于“我這里”是多大的范圍，根據(jù)公里級(jí)的IP地理信息應(yīng)該可以推測(cè)出來(lái)。這不是本文的重點(diǎn)。

本題在討論過(guò)程中，得到了色影無(wú)忌眾多網(wǎng)友的見(jiàn)解，@siliconworm, @TakeThat05, @rks, @pipedream，特此致謝。