博客園排名預測
前言
之前寫過一篇繪製博客園積分與排名趨勢圖的文章——《查看博客園積分與排名趨勢圖的工具 》，使用那篇文章介紹的工具，可以通過趨勢圖直觀的看出排名前進的走勢。但是如果想看看自己積分達到多少才能進入前多少名次，就無能為力了。如果我們能夠根據歷史數據，擬合出一條預測曲線，然後根據這條曲線就可以預測多少積分進多少排名啦！想想就很激動吶~
積分-排名曲線
在開始擬合數據之前，我們先看下現在的趨勢圖：
它的橫軸是時間軸，兩條縱軸分別是積分與排名。雖然積分總是增長、排名大體是下降趨勢，但是這個趨勢我試了幾種方法都不太好擬合，感覺沒什麼規律，特別受發表文章時機的影響：如果每天發一篇，這個曲線肯定增長的快；如果一年也不發，增長的肯定慢。我們的預測其實和時間沒有什麼關係，主要自變量是積分，因變量是排名。那麼是否可以做一張圖，橫軸是積分，縱軸是排名呢？於是就有了下面的改進：
這下清爽多了，可以看到，雖然發表文章時積分會快速增長從而在橫軸留下一些空白，但是總的來看積分與排名相對趨勢是維持在一條曲線上的。想要繪製這樣一條曲線，gnuplot 腳本改動並不大：
1 #! /usr/bin/gnuplot 2 set terminal png size 1080,720 #建立空白圖片 3 set title usr.": score (".y1range.") rank (".y2range.")" #註明曲線圖標題 4 set output "fit.png" #設置文件名 5 set key left 6 set grid 7 8 set xlabel "score" 9 set ylabel "rank" 10 11 #plot "score.txt" using 1:2 with lp pt 13 title "score" axis x1y1, "score.txt" using 1:3 with lp pt 3 title "rank" axis x1y2 12 plot "score.txt" using 2:3 with lp pt 13 title "score-rank" 13 14 quit #退出軟件
就是將 line 11 換成 line 12，再去掉一些冗餘代碼就可以了。
數據擬合
有了歷史數據和正確的映射關係，就可以進行數據擬合了。數據擬合最重要的是找到擬合函數，第一眼看到上面那條曲線我想到的就是二次函數，可以用拋物線的一段來進行擬合。此外隨着積分的無限增大，排名最終會無限貼進於 1，這個特性讓我想到了很多無限貼進某個值的函數，例如倒數函數無限貼進於 0，對數函數的反函數無限貼近於 0，還有反正切函數也是如此。最近剛剛高考結束，相信大家對高中數學已經忘了個七七八八，還好我這裡準備了幾張圖：
從左到右依次是倒數、對數、反正切(atan)，對於後兩者，還需要加工一下才能滿足我們的要求，具體羅列如下：
f1(x)=a*x^2+b*x+c f2(x)=f/x+g f3(x)=j*atan(x)+k f4(x)=m*log(x)+n
分別就是它們四個啦，為了之後區分參數方便，使用了不同的參數名稱。有的人可能會問了，這個對數是無限趨近於無窮大的，怎麼能和上面的曲線擬合在一起呢？反正切也存在類似的問題，先別著急，讓我們看下擬合結果就知道了。
嘿，這樣居然也行，四條曲線竟然都能擬合上。gnuplot 腳本改動並不大：
1 #! /usr/bin/gnuplot 2 set terminal png size 1080,720 #建立空白圖片 3 set title usr.": score (".y1range.") rank (".y2range.")" #註明曲線圖標題 4 set output "fit.png" #設置文件名 5 set key left 6 set grid 7 8 set xlabel "score" 9 set ylabel "rank" 10 11 y1(x)=a*x**2+b*x+c 12 fit y1(x) "score.txt" using 2:3 via a,b,c 13 14 y2(x)=f/x+g 15 fit y2(x) "score.txt" using 2:3 via f,g 16 17 y3(x)=j*atan(x)+k 18 fit y3(x) "score.txt" using 2:3 via j,k 19 20 y4(x)=m*log(x)+n 21 fit y4(x) "score.txt" using 2:3 via m,n 22 23 plot "score.txt" using 2:3 with lp pt 13 title "score-rank", \ 24 y1(x) with l lw 2 lt 2 title "f(x)=ax^2+bx+c", \ 25 y2(x) with l lw 3 lt 3 title "f(x)=a/x+b", \ 26 y3(x) with l lw 1 lt 4 title "f(x)=a*atan(x)+b", \ 27 y4(x) with l lw 2 lt 5 title "f(x)=a*log(x)+b" 28 29 quit #退出軟件
line 11-21 定義了各個擬合函數，line 24-27 增加了擬合曲線的繪製。如果能將擬合後的函數參數標識出來，就更好了，其實也不難，因為 a/b/c/f/g/j/k/m/n 這些參數在 gnuplot 腳本中就可以直接訪問，只需要在圖例顯示處增加一些代碼就可以了：
plot "score.txt" using 2:3 with lp pt 13 title "score-rank", \ y1(x) with l lw 6 lt 2 title sprintf("f1(x)=%.8fx^2%+fx%+.0f",a,b,c), \ y2(x) with l lw 5 lt 3 title sprintf("f2(x)=%.2f/x%+.0f",f,g), \ y3(x) with l lw 2 lt 4 title sprintf("f3(x)=%.2fatan(x)%+.0f",j,k), \ y4(x) with l lw 2 lt 5 title sprintf("f4(x)=%.2flog(x)%+.0f",m,n)
最終效果如下：
看到對數函數的參數，可以解答之前的疑問了，m=-39186.57 (<0) 使得整個對數函數走勢是從高到低的，從而符合歷史數據的趨勢，反正切函數也是如此。另外從圖中還可以看到另一個有意思的現象，就是 f2 倒數函數和 f3 反正切函數完全重合！為了突出這一點，我特意加粗了 f2，變細了 f3，這樣就不會像之前那樣看不清楚了。觀察兩個函數的參數，乘法係數 f 和 j 非常接近，加法係數 g 和 k 是非常不一樣的，最終卻異曲同工走在了一條路上，真是不可思議。某博士同學告訴我反正切就是倒數，所以會重合，原理我沒聽懂，不過既然這兩個曲線一樣，那就保留一樣好了，這裡我選擇了倒數函數，相對直觀一些。
數據預測
有了三個擬合函數，就可以對數據進行預測了，一開始雄心勃勃，打算預測一下自己 40 W 分時的排名 (有點不自量力哈)，預測值通過 label 形式輸出在圖形上，就像這樣：
結果相去甚遠，首先恭喜二次函數 (f1) 得到了 2800 W 名的好成績~ 話說我剛加入博客園也才 14 W 名，簡直離譜；倒數函數得到了 4500 名的結果，比現在 (5900) 也好不了多少，這 38 W 分看來白漲了； 對數函數最牛，直接干成負數了，祝賀我進入了另一個維度~ 看來還是我好高騖遠了，手裡拿着 2 W 分的數據，卻想預測 40 W 分的排名，難為 gnuplot 了。調整了一下目標，先算個 4 W 分的吧：
現在稍微靠譜一些了，總的來說，二次函數偏差太大，倒數函數偏大，對數函數偏小。其實想想也能明白，二次函數拋物線嘛，遲早要轉回來的，偏大也可以理解。最後上 gnuplot 腳本：
1 #! /usr/bin/gnuplot 2 set terminal png size 1080,720 #建立空白圖片 3 set title usr.": score (".y1range.") rank (".y2range.")" #註明曲線圖標題 4 set output "fit.png" #設置文件名 5 set key left 6 set grid 7 8 set xlabel "score" 9 set ylabel "rank" 10 11 y1(x)=a*x**2+b*x+c 12 fit y1(x) "score.txt" using 2:3 via a,b,c 13 14 y2(x)=f/x+g 15 fit y2(x) "score.txt" using 2:3 via f,g 16 17 y3(x)=m*log(x)+n 18 fit y3(x) "score.txt" using 2:3 via m,n 19 20 xval=40000 21 y1val=a*xval**2+b*xval+c 22 y2val=f/xval+g 23 y3val=m*log(xval)+n 24 25 set label 1 sprintf("f1(%.0f)=%.0f",xval,y1val) at graph 0.6,0.7 left 26 set label 2 sprintf("f2(%.0f)=%.0f",xval,y2val) at graph 0.6,0.65 left 27 set label 3 sprintf("f3(%.0f)=%.0f",xval,y3val) at graph 0.6,0.6 left 28 29 plot "score.txt" using 2:3 with lp pt 13 title "score-rank", \ 30 y1(x) with l lw 4 lt 2 title sprintf("f1(x)=%.8fx^2%+fx%+.0f",a,b,c), \ 31 y2(x) with l lw 3 lt 3 title sprintf("f2(x)=%.2f/x%+.0f",f,g), \ 32 y3(x) with l lw 2 lt rgb "red" title sprintf("f3(x)=%.2flog(x)%+.0f",m,n) 33 34 35 quit #退出軟件
主要增加的內容就是中間幾行啦，line 20-23 通過現有參數計算 x 為 40000 時的各函數 y 值，line 25-27 將計算好的預測值顯示在 label 中。其實函數已經定義好了，如果能直接通過 f1 (40000) / f2 (40000) / f3 (40000) 得到結果就更好了，但是沒有在 gnuplot 手冊中找到這種語法，不得己自己再寫一遍，有懂行的同學不吝賜教哈。最後因為我們的預測值都是整數，所以打印出來的數據也沒有保留小數位，通過 sprintf 自動四捨五入了。
繪製預測曲線
上面的代碼可以預測某個點的數據，但是還是有點呆板，需要手動指定預測值，如果將預測值設置為當前分數的兩倍，就能自動預測啦。將得到的預測值寫入一個數據文件，隨着時間積累，形成一條預測曲線繪製出來，再和實際數據做對比，預測效果豈不一目了然？
輸出預測值
將 gnuplot 腳本中計算得到的預測值寫入一個文件，這個事情看起來簡單做起來難，難就難在我找了半天，沒有找到可以從腳本直接輸出信息到 console 或重定向到文件的方法。echo 這種命令在 gnuplot 腳本中是不存在的，於是這裡繞了一個大圈——在腳本執行完成後，通過分拆 fit.log 中的擬合日誌提取函數的各個參數 (a/b/c/f/g/m/n)，再構建函數計算預測值，最後寫入數據文件——哪位高手如果知道如何在 gnuplot 腳本中直接輸出信息的話，不吝賜教哈，就可以把這個大彎路省掉了。
現在來看這個蹩腳的彎路，在開始參數提取前，先熟悉一下 gnuplot 的擬合日誌：
*******************************************************************************
Mon Jul 5 14:22:31 2021
FIT: data read from "score.txt" using 2:3
format = x:z
#datapoints = 365
residuals are weighted equally (unit weight)
function used for fitting: y1(x)
y1(x)=a*x**2+b*x+c
fitted parameters initialized with current variable values
iter chisq delta/lim lambda a b c
0 1.2882895936e+19 0.00e+00 1.08e+08 1.000000e+00 1.000000e+00 1.000000e+00
11 6.0109804999e+08 -9.02e-01 1.08e-03 1.991252e-04 -8.363164e+00 1.465348e+05
After 11 iterations the fit converged.
final sum of squares of residuals : 6.01098e+08
rel. change during last iteration : -9.01929e-06
degrees of freedom (FIT_NDF) : 362
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 1288.6
variance of residuals (reduced chisquare) = WSSR/ndf : 1.66049e+06
Final set of parameters Asymptotic Standard Error
======================= ==========================
a = 0.000199125 +/- 3.57e-06 (1.793%)
b = -8.36316 +/- 0.08565 (1.024%)
c = 146535 +/- 456.8 (0.3117%)
correlation matrix of the fit parameters:
a b c
a 1.000
b -0.986 1.000
c 0.921 -0.969 1.000
*******************************************************************************
Mon Jul 5 14:22:31 2021
FIT: data read from "score.txt" using 2:3
format = x:z
#datapoints = 365
residuals are weighted equally (unit weight)
function used for fitting: y2(x)
y2(x)=f/x+g
fitted parameters initialized with current variable values
iter chisq delta/lim lambda f g
0 2.5365572521e+12 0.00e+00 7.07e-01 1.000000e+00 1.000000e+00
7 2.5241195880e+09 -4.95e-08 7.07e-08 3.478261e+08 4.432964e+04
After 7 iterations the fit converged.
final sum of squares of residuals : 2.52412e+09
rel. change during last iteration : -4.94761e-13
degrees of freedom (FIT_NDF) : 363
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 2636.95
variance of residuals (reduced chisquare) = WSSR/ndf : 6.9535e+06
Final set of parameters Asymptotic Standard Error
======================= ==========================
f = 3.47826e+08 +/- 2.776e+06 (0.7982%)
g = 44329.6 +/- 327.3 (0.7383%)
correlation matrix of the fit parameters:
f g
f 1.000
g -0.907 1.000
*******************************************************************************
Mon Jul 5 14:22:31 2021
FIT: data read from "score.txt" using 2:3
format = x:z
#datapoints = 365
residuals are weighted equally (unit weight)
function used for fitting: y3(x)
y3(x)=m*log(x)+n
fitted parameters initialized with current variable values
iter chisq delta/lim lambda m n
0 2.5360128666e+12 0.00e+00 6.58e+00 1.000000e+00 1.000000e+00
5 2.4184679129e+08 -5.46e-07 6.58e-05 -3.918657e+04 4.438066e+05
After 5 iterations the fit converged.
final sum of squares of residuals : 2.41847e+08
rel. change during last iteration : -5.46492e-12
degrees of freedom (FIT_NDF) : 363
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 816.238
variance of residuals (reduced chisquare) = WSSR/ndf : 666245
Final set of parameters Asymptotic Standard Error
======================= ==========================
m = -39186.6 +/- 95.82 (0.2445%)
n = 443807 +/- 886.9 (0.1998%)
correlation matrix of the fit parameters:
m n
m 1.000
n -0.999 1.000
它將擬合的迭代、最終的參數、誤差率清楚的羅列了出來，這裡我們需要提取的只是各個參數最終的值。好在這些行都比較有特點，基本遵循 “參數名 = 數據 ***” 的格式，於是可以用 grep 先過濾一把：
sed -n '/[abcfgmn] *=.*/p' fit.log
經過這樣處理後，就只剩下這樣的內容了：
a = 0.000199125 +/- 3.57e-06 (1.793%) b = -8.36316 +/- 0.08565 (1.024%) c = 146535 +/- 456.8 (0.3117%) f = 3.47826e+08 +/- 2.776e+06 (0.7982%) g = 44329.6 +/- 327.3 (0.7383%) m = -39186.6 +/- 95.82 (0.2445%) n = 443807 +/- 886.9 (0.1998%)
再使用 awk 提取我們關心的前三列：
sed -n '/[abcfgmn] *=.*/p' fit.log | awk '{print $1,$2,$3+0}'
注意第三列使用 “$3+0” 的 trick 來保證提取的是浮點數據：
a = 0.000199125 b = -8.36316 c = 146535 f = 347826000 g = 44329.6 m = -39186.6 n = 443807
最後將它們傳遞給 eval 實現「求值」，即當前 shell 中自動獲得了 a-n 7 個變量和它們的初值：
eval $(sed -n '/[abcfgmn] *=.*/p' fit.log | awk '{print $1,$2,$3+0}' | sed 's/ //g')
注意 shell 中的賦值語句是不能有空格的，所以需要使用 sed 過濾一下空格。有了函數的各個參數，就可以利用函數計算預測值了，這裡我使用了 awk，因為它對浮點數運算支持的比較好：
awk -v a=$a -v b=$b -v c=$c -v f=$f -v g=$g -v m=$m -v n=$n -v xval=$xval 'BEGIN { print "y1="int(a*xval*xval+b*xval+c+0.5); print "y2="int(f/xval+g+0.5); print "y3="int(m*log(xval)+n+0.5) }'
這裡首先利用 awk 的 -v 選項將 shell 腳本中的變量傳遞到 awk 中，然後在 awk 中根據三個函數分別計算了三個預測值。xval 是事先定義好的，一般就是當前 x 值的 2 倍。上面的腳本輸出如下：
y1=130609 y2=53025 y3=28561
這裡的四捨五入使用了 +0.5 的笨辦法，最終結果和 gnuplot 計算的完全一致。有了這個就可以繼續利用 eval 搞事情了，下面上完整的提取、計算、寫入腳本：
1 #! /bin/sh 2 usr=$(cat user.txt) 3 firstline=$(cat ./score.txt | head -1) 4 y1min=$(echo $firstline | awk '{ print $2 }') 5 y2max=$(echo $firstline | awk '{ print $3 }') 6 lastline=$(cat ./score.txt | tail -1) 7 y1max=$(echo $lastline | awk '{ print $2 }') 8 y2min=$(echo $lastline | awk '{ print $3 }') 9 echo "y1 range [$y1min,$y1max], y2 range [$y2max,$y2min]" 10 gnuplot -e "usr='$usr'" -e "y1range='$(($y1max-$y1min))'" -e "y2range='$(($y2max-$y2min))'" ./draw.plt 11 12 # clear log to get generated values 13 rm fit.log 2>/dev/null 14 gnuplot -e "usr='$usr'" -e "y1max='$y1max'" -e "y2min='$y2min'" ./fit.plt 15 16 if [ $# -gt 0 ]; then 17 # don't know how to print out parameter (a,b,c,f,g,m,n) from gnuplot script, so here extract them from fit.log (that's why we clear fit.log before calling fit.plt) 18 eval $(sed -n '/[abcfgmn] *=.*/p' fit.log | awk '{print $1,$2,$3+0}' | sed 's/ //g') 19 # after that line, a/b/c/f/g/m/n takes effect, now calculate predicating values 20 21 # predicate x*2, round result to integer 22 xval=$(($y1max*2)) 23 eval $(awk -v a=$a -v b=$b -v c=$c -v f=$f -v g=$g -v m=$m -v n=$n -v xval=$xval 'BEGIN { print "y1="int(a*xval*xval+b*xval+c+0.5); print "y2="int(f/xval+g+0.5); print "y3="int(m*log(xval)+n+0.5) }') 24 25 # dump results to data files 26 echo "$xval $y1" >> predicate_binomial.data 27 echo "$xval $y2" >> predicate_reciprocal.data 28 echo "$xval $y3" >> predicate_logarithm.data 29 fi 30 31 # for centos 32 type eog > /dev/null 2>&1 33 if [ $? -eq 0 ]; then 34 eog draw.png & 35 eog fit.png & 36 exit 0 37 fi 38 39 # for mac 40 type open > /dev/null 2>&1 41 if [ $? -eq 0 ]; then 42 open draw.png & 43 open fit.png & 44 exit 0 45 fi 46 47 # for windows msys2 48 type mspaint > /dev/null 2>&1 49 if [ $? -eq 0 ]; then 50 mspaint draw.png & 51 mspaint fit.png & 52 exit 0 53 fi 54 55 exit 1
重點就是 line 16-29 啦，這段代碼只有在腳本參數個數大於 0 時才生效，也就是說你平時做一個 ./plot.sh 這種動作查看趨勢圖時，是不會修改數據的，只有當定時任務每日執行 score.sh 才會調用帶參數的 plot.sh 去更新預測值。關於 score.sh 的內容，可以參數我之前寫的那篇文章。
預測值經過計算並提取到 shell 腳本後，分別存儲在了三個 data 文件中，文件名說明了他們使用的擬合函數。這裡也可以將多個擬合函數的預測值放在一個文件，畢竟他們的 x 軸數據都是一樣的嘛，沒有這樣做的原因主要是考慮到後期可能加入新的擬合函數來進行預測，獨立存儲的話互不干擾，加入刪除都比較方便，利於擴展。
歷史補算
以現在 20000 分得到的預測值是在 40000 分左右，想要驗證預測值準不準，還需要我的博客漲到 40000 分才能拉在一起做對比，以我現在這種速度，不知道要等到猴年馬月，然而我是這種沒有耐心的人嗎？確實是！於是我打算將 5000-20000 這段歷史數據的預測值也給它補全了，這樣生成的預測值就涵蓋了 10000-40000 的範圍，正好能和我 10000-20000 的區間形成交集，從而對比着看。
說起來複雜，做起來簡單，我只需將歷史數據一分為二，一部分是 4000-5000 的極小規模的歷史數據；一部分是 5000-20000 的補算歷史數據。使用 plot.sh 作用於第一部分數據，生成預測值，然後從第二部分數據頭部取出一條記錄添加到第一部分數據末尾，再調用 plot.sh 生成一條預測數據……周而復始，直到第二部分數據消耗完畢。下面就是用於補算的腳本：
1 #! /bin/sh 2 # data should be parted into score.txt (with some basic historical data) 3 # and more.txt (with more recent data), 4 # and then we will repair predicating data by moving more.txt data 5 # into score.txt line by line and calling plot.sh without call png starter... 6 7 if [ ! -f "more.txt" ]; then 8 echo "you should split score.txt to score.txt & more.txt first before run this scripts..." 9 exit 1 10 fi 11 12 while read line 13 do 14 echo "repair line $line" 15 echo "$line" >> score.txt 16 ./plot.sh "update_predicating_data" 17 done < more.txt 18 19 echo "repair done, now submit predicat_*.data !" 20 rm more.txt
當然了，在運行這個腳本之前，你需要將 score.txt 一切為二，前一部分還叫 score.txt，第二部分需要命名為 more.txt。當腳本運行完畢後，more.txt 中的數據自動合入 score.txt，同時產生三個包含歷史預測值的 predicate_xxx.data 文件。而且為了防止在補算歷史過程中自動打開大量的趨勢圖文件，需要在 plot.sh 中做如下修改：
1 if [ $# -gt 0 ]; then 2 # don't know how to print out parameter (a,b,c,f,g,m,n) from gnuplot script, so here extract them from fit.log (that's why we clear fit.log before calling fit.plt) 3 eval $(sed -n '/[abcfgmn] *=.*/p' fit.log | awk '{print $1,$2,$3+0}' | sed 's/ //g') 4 # after that line, a/b/c/f/g/m/n takes effect, now calculate predicating values 5 6 # predicate x*2, round result to integer 7 xval=$(($y1max*2)) 8 eval $(awk -v a=$a -v b=$b -v c=$c -v f=$f -v g=$g -v m=$m -v n=$n -v xval=$xval 'BEGIN { print "y1="int(a*xval*xval+b*xval+c+0.5); print "y2="int(f/xval+g+0.5); print "y3="int(m*log(xval)+n+0.5) }') 9 10 # dump results to data files 11 echo "$xval $y1" >> predicate_binomial.data 12 echo "$xval $y2" >> predicate_reciprocal.data 13 echo "$xval $y3" >> predicate_logarithm.data 14 else 15 # for centos 16 type eog > /dev/null 2>&1 17 if [ $? -eq 0 ]; then 18 eog draw.png & 19 eog fit.png & 20 exit 0 21 fi 22 23 # for mac 24 type open > /dev/null 2>&1 25 if [ $? -eq 0 ]; then 26 open draw.png & 27 open fit.png & 28 exit 0 29 fi 30 31 # for windows msys2 32 type mspaint > /dev/null 2>&1 33 if [ $? -eq 0 ]; then 34 mspaint draw.png & 35 mspaint fit.png & 36 exit 0 37 fi 38 fi 39 40 exit 1
line 14 添加了一行 else，表示只有在腳本參數為 0 時才啟動圖片自動打開功能。
繪製預測線
前面鋪墊了這麼多，終於可以把預測值繪製出來一睹芳容了：
先撇開預測曲線的風騷走位，重點關注一下 10000-20000 這個區間，可以看到點劃線的真實數據和三條預測曲線相差還是蠻大的，只有對數函數在一開始非常貼近真實值，後來也出現偏離。不過總的來說對數最准、倒數整體偏高、二次函數就是來搞笑的——上下橫跳沒準頭。有的人可能奇怪了，這個預測值為什麼和擬合曲線差距這麼大？原因是預測曲線的每個點的參數都不一樣，由之前小一半的歷史數據擬合計算得到的，所以不能完美重合擬合函數，可以將預測曲線理解成是一堆擬合函數的末位點集合形成的軌跡 (稍費腦，理解不了就不用理解了)。最後再來欣賞一下全圖，無意間拉大橫軸範圍後，二次函數的拋物線本性果然暴露了，對數和倒數表現倒還可以。過足了預測癮後，還是讓我們限制一下橫軸範圍，不然真實數據實在是看不清了：
同時調整的還有圖例顯示方式和預測曲線的顏色，後者可以讓你一眼看出擬合函數與預測曲線關係。下面是最終的 gnuplot 腳本：
1 #! /usr/bin/gnuplot 2 set terminal png size 1080,720 #建立空白圖片 3 set title usr.": score (".y1max.") rank (".y2min.")" #註明曲線圖標題 4 set output "fit.png" #設置文件名 5 set key left reverse Left spacing 1.2 6 set grid 7 8 set xlabel "score" 9 set ylabel "rank" 10 # to prevent predicating value pollute our x-axis 11 set xrange [y1min-100:y1max+100] 12 # no fit log in console (fit.log only) 13 set fit quiet 14 15 y1(x)=a*x**2+b*x+c 16 fit y1(x) "score.txt" using 2:3 via a,b,c 17 18 y2(x)=f/x+g 19 fit y2(x) "score.txt" using 2:3 via f,g 20 21 y3(x)=m*log(x)+n 22 fit y3(x) "score.txt" using 2:3 via m,n 23 24 xval=y1max*2 25 y1val=a*xval**2+b*xval+c 26 y2val=f/xval+g 27 y3val=m*log(xval)+n 28 29 set label 1 sprintf("f1(%.0f)=%.0f",xval,y1val) at graph 0.6,0.7 left 30 set label 2 sprintf("f2(%.0f)=%.0f",xval,y2val) at graph 0.6,0.65 left 31 set label 3 sprintf("f3(%.0f)=%.0f",xval,y3val) at graph 0.6,0.6 left 32 33 plot "score.txt" using 2:3 with lp pt 13 title "score-rank", \ 34 y1(x) with l lw 4 lt 2 title sprintf("f1(x)=%.8fx^2%+fx%+.0f",a,b,c), \ 35 y2(x) with l lw 3 lt 3 title sprintf("f2(x)=%.2f/x%+.0f",f,g), \ 36 y3(x) with l lw 2 lt rgb "red" title sprintf("f3(x)=%.2flog(x)%+.0f",m,n), \ 37 "predicate_binomial.data" using 1:2 with lp pt 12 lt 2 title "f1-pred", \ 38 "predicate_reciprocal.data" using 1:2 with lp pt 11 lt 3 title "f2-pred", \ 39 "predicate_logarithm.data" using 1:2 with lp pt 10 lt rgb "red" title "f3-pred" 40 41 quit #退出軟件
主要改動點在 line 37-39，使用之前補算的預測值繪圖。
結論
最後比較靠譜的預測值就是取對數函數與倒數函數預測值靠近對數函數 1/3 的位置，即 f(x)=(f3(x)-f2(x))/3+f2(x)=2/3*f2(x)+1/3*f3(x)。
這個就是我胡謅的啦~ 真正的結論就是自己的數值分析沒學好，需要回爐重造，連一個好的擬合函數也找不到，此刻終於在畢業 N 多年後認識到高等數據和數值分析的重要性…
後記
文中代碼可以從這裡下載：
//github.com/goodpaperman/cnblogs
這個代碼庫可以直接 fork 喲~ 在寫文章過程中發現了一個很有意思的數學教學網站，可以在 web 端直接繪製你輸入的表達試，推薦一波：//www.shuxuele.com/index.html
參考
[1]. gnuplot圖例legend設置
[2]. awk將字符串轉為數字的方法
[3]. 在命令行中使用gnuplot快速查看數據
[4]. Gnuplot重定向fit輸出
[5]. gnuplot常用技巧
[6]. 在gnuplot中，繪製一些分段函數
[7]. gnuplot使用手冊
[9]. AWK 打印匹配內容之後的指定行
