這是一篇正經(jīng)的數(shù)據(jù)分析案例。
去年12月初，在經(jīng)過四年多的積累后，編程教室微信公眾號的關(guān)注人數(shù)突破10萬人。（可回顧 最開始我也沒有想過會(huì)有這么一天…）
10萬人只是另一個(gè)開始，讓我感到責(zé)任更大了。如果不寫出更多更好質(zhì)量的文章和教程，也對不起大家的關(guān)注啊。人數(shù)不是目的，內(nèi)容才是王道。
但是嘛，偶爾也會(huì) yy 一下，什么時(shí)候我們的關(guān)注數(shù)能到達(dá)更高的量級，比如，100萬？
既然 Python 可以用來做數(shù)據(jù)分析，何不根據(jù)我們公眾號現(xiàn)有的用戶增長數(shù)據(jù)來分析一下，什么時(shí)候可以迎來第100萬個(gè)關(guān)注者？
說干就干！（不想看過程的直接拉到最后看結(jié)論）
微信后臺(tái)可以導(dǎo)出用戶增長數(shù)據(jù)的 excel 表格。數(shù)據(jù)從2013年7月開始，每次導(dǎo)出時(shí)間間隔最多一年。我們編程教室的賬號是2013年6月份創(chuàng)建的，雖然差了一點(diǎn)點(diǎn)，但也足夠了。
把幾年的數(shù)據(jù)合并一下，我們這次只關(guān)注其中的“累積關(guān)注人數(shù)”和“時(shí)間”。通過 matplotlib 把增長曲線繪制出來
顯然這不是一個(gè)簡單的勻速增長曲線，而是加速增長！這讓我甚感欣慰。
核心代碼
# x_data 時(shí)間列表
# y_data 用戶數(shù)列表
plt.plot(x_data, y_data,'g')
plt.show()
那么從數(shù)學(xué)上來看，有沒有能夠較好擬合這個(gè)增長曲線的模型呢？我們來嘗試幾個(gè)最常用的擬合曲線，看看效果。
多項(xiàng)式擬合即用形如
的函數(shù)曲線來擬合現(xiàn)有的數(shù)據(jù)。比如三次多項(xiàng)式擬合就是對公式
中的4個(gè)系數(shù)求解，使得函數(shù)曲線與數(shù)據(jù)“最接近”。
至于怎樣才算是“最接近”？直觀來考慮，就是擬合曲線和實(shí)際曲線上對應(yīng)點(diǎn)的距離最短，即絕對值最小。以我們的例子來說，就是擬合函數(shù)算出的每天總關(guān)注人數(shù)和當(dāng)天實(shí)際總關(guān)注人數(shù)的差，我們要讓這個(gè)差值的總和最小。
但因?yàn)榻^對值之和不容易處理，所以通常我們選擇差值的平分和來替代。這就是“最小二乘法”。
更數(shù)學(xué)化的表述就是，我們要找出擬合曲線中的一組參數(shù) c，使得模型與實(shí)際值上每一點(diǎn)的殘差 ek 的平方和最小。
我們繪制了從1次多項(xiàng)式（線性函數(shù)）到9次多項(xiàng)式的擬合曲線：
因?yàn)槲覀兊哪康氖且A(yù)測之后的趨勢，所以選擇的擬合天數(shù)要大于實(shí)際數(shù)據(jù)的天數(shù)。
從圖上就能比較直觀地就看出，1次、2次等低階曲線不能很好地貼合原數(shù)據(jù)，3~8次效果都還不錯(cuò)，而9次曲線在不久之后就會(huì)因?yàn)檫^擬合而產(chǎn)生不合理的波動(dòng)。
對于多項(xiàng)式擬合，numpy 提供了現(xiàn)成的 polyfit 和 poly1d 函數(shù)供調(diào)用。
核心代碼
# x_np 實(shí)際數(shù)據(jù)，時(shí)間
# y_np 實(shí)際數(shù)據(jù)，用戶數(shù)
# x_fit 擬合數(shù)據(jù)，時(shí)間
coeff = np.polyfit(x_np, y_np, k)
poly = np.poly1d(coeff)
y_fit = poly(x_fit)
plt.plot(x_fit, y_fit)
指數(shù)函數(shù)是重要的基本初等函數(shù)之一，這里我們通過確定以 e 為底的函數(shù)
中3個(gè)參數(shù) a、b、c 來進(jìn)行擬合。
看起來擬合效果還不錯(cuò)。
numpy 沒有提供直接的指數(shù)擬合函數(shù)，但我們可以通過 scipy 庫里的 scipy.optimize.leastsq 實(shí)現(xiàn)最小二乘法。
核心代碼
def func(x, p):
a,b,c = p
return a * np.exp(b * x) + c
# 殘差函數(shù)
def residuals(p, y, x):
return y - func(x, p)
pe = [1, 0.0001, 1] # 初始預(yù)測值
plsq = leastsq(residuals, pe, args=(y_np, x_np))
y_fit = func(x_fit, plsq[0])
plt.plot(x_fit, y_fit)
冪函數(shù)和指數(shù)函數(shù)有點(diǎn)類似，只不過我們使用的函數(shù)是
同樣也是3個(gè)參數(shù)。
擬合的效果與前面的指數(shù)函數(shù)有點(diǎn)相似。代碼中，我們也只要在剛才的基礎(chǔ)上，修改一下 func 函數(shù)即可。
核心代碼
def func(x, p):
a,b,c = p
return a * x ** b + c
以上幾種方法雖然看起來都不錯(cuò)，但結(jié)果畢竟有不小的差異，究竟哪一個(gè)更“科學(xué)”一點(diǎn)呢？
我們通過幾個(gè)評價(jià)指標(biāo)來衡量一下：
均方根誤差（RMSE）：真實(shí)值和預(yù)測值之差的平方和。這其實(shí)就是我們擬合時(shí)的判斷基礎(chǔ)啊。只不過加上了根號，使得結(jié)果的量綱更加合理（否則就是均方誤差MSE）。
平均絕對誤差（MAE）：和 MSE 的區(qū)別就在于直接使用真實(shí)值和預(yù)測值之差的絕對值作為衡量標(biāo)準(zhǔn)。
R平方（R2） ：因?yàn)?MSE 結(jié)果的大小取決于不同數(shù)據(jù)的本身數(shù)值大小，并不統(tǒng)一。R2 則是在此基礎(chǔ)上，將其轉(zhuǎn)換至 0~1 之間，以便于評價(jià)。
以上指標(biāo)，sklearn 庫均在 metrics 中提供了方法。
核心代碼
# ploy 擬合函數(shù)
rmse = sqrt(metrics.mean_squared_error(y_np, poly(x_np)))
mae = metrics.mean_absolute_error(y_np, poly(x_np))
r2 = metrics.r2_score(y_np, poly(x_np))
當(dāng)然，這些指標(biāo)都是基于擬合函數(shù)與已有數(shù)據(jù)的判斷，對于未來的預(yù)測，誰也說不準(zhǔn)，只能是“僅供參考”。畢竟如果可以預(yù)知未來，那我大概早就 all in 比特幣了。
| 函數(shù) | 100萬用戶 | RMSE | MAE | R2 |
1次 | 2063/4/15 | 12132 | 9388 | 0.846 |
2次 | 2026/5/18 | 4802 | 3377 | 0.974 |
3次 | 2022/6/25 | 2939 | 1765 | 0.990 |
4次 | 2021/3/12 | 2637 | 1957 | 0.992 |
5次 | 2019/12/12 | 1901 | 1519 | 0.996 |
6次 | 2019/5/16 | 1143 | 748 | 0.999 |
7次 | 2019/3/6 | 994 | 507 | 0.999 |
8次 | 2019/3/11 | 994 | 510 | 0.999 |
9次 | — | 975 | 531 | 0.999 |
指數(shù) | 2020/12/16 | 2682 | 1722 | 0.992 |
冪函數(shù) | 2022/12/28 | 3462 | 2440 | 0.987 |
綜合結(jié)果來看，編程教室的百萬用戶很可能在2019~2022年之間到來。對3~8次多項(xiàng)式、指數(shù)函數(shù)、冪函數(shù)的預(yù)測結(jié)果做個(gè)簡單的平均，那么這一天就是：
2020年5月27日
只需要 811 天，想想還有點(diǎn)小激動(dòng)呢。
忽然，我想到了那個(gè)詭異的9次函數(shù)，說來也不是不可能哦：當(dāng)人數(shù)過了40萬，因?yàn)槟硞€(gè)不小心被封了號，一切歸零。這也不是什么新鮮事兒。
所以，我還是老老實(shí)實(shí)寫教程吧。猥瑣發(fā)育，別浪！
聯(lián)系客服
