前言
之前写过一篇绘制博客园积分与排名趋势图的文章——《查看博客园积分与排名趋势图的工具 》，使用那篇文章介绍的工具，可以通过趋势图直观的看出排名前进的走势。但是如果想看看自己积分达到多少才能进入前多少名次，就无能为力了。如果我们能够根据历史数据，拟合出一条预测曲线，然后根据这条曲线就可以预测多少积分进多少排名啦！想想就很激动呐~
积分-排名曲线
在开始拟合数据之前，我们先看下现在的趋势图：
它的横轴是时间轴，两条纵轴分别是积分与排名。虽然积分总是增长、排名大体是下降趋势，但是这个趋势我试了几种方法都不太好拟合，感觉没什么规律，特别受发表文章时机的影响：如果每天发一篇，这个曲线肯定增长的快；如果一年也不发，增长的肯定慢。我们的预测其实和时间没有什么关系，主要自变量是积分，因变量是排名。那么是否可以做一张图，横轴是积分，纵轴是排名呢？于是就有了下面的改进：
这下清爽多了，可以看到，虽然发表文章时积分会快速增长从而在横轴留下一些空白，但是总的来看积分与排名相对趋势是维持在一条曲线上的。想要绘制这样一条曲线，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 多年后认识到高等数据和数值分析的重要性…
后记
文中代码可以从这里下载：
https://github.com/goodpaperman/cnblogs
这个代码库可以直接 fork 哟~ 在写文章过程中发现了一个很有意思的数学教学网站，可以在 web 端直接绘制你输入的表达试，推荐一波：https://www.shuxuele.com/index.html
参考
[1]. gnuplot图例legend设置
[2]. awk将字符串转为数字的方法
[3]. 在命令行中使用gnuplot快速查看数据
[4]. Gnuplot重定向fit输出
[5]. gnuplot常用技巧
[6]. 在gnuplot中，绘制一些分段函数
[7]. gnuplot使用手册
[9]. AWK 打印匹配内容之后的指定行
[10]. UWP 获取博客园积分，并以图表形式呈现变化趋势