用上市公司2019年报净利润数据验证本福特定律

在疫情期间，听到一个消息：有国外研究者用“本福特”定律验证了中国公布的疫情信息，显示数据是真实可信的。这是第一次听到这个名字，感觉很神奇，它竟然能够验证数据有无作假。
上网查了一下，原来所谓本福特定律是指在自然状态下形成的一组数据，数字的首位分布应当符合一定的规律。具体来说从1到9的9个阿拉伯数字i，出现在首位的概率应当是：P(i)=log10((i+1)/i)。
我于是很感兴趣，想用股票数据来验证一下这个本福特定律。
我们采用tushare接口获取2019年年报（第4季度）数据，取其中的净利润数据，然后我们只考虑净利润为正的情况。

# 验证本福特定律
import tushare as ts 
import math
import matplotlib.pyplot as plt
import pandas as pd
from functools import reduce
from pylab import *
# 这一句让pyplot支持中文显示
mpl.rcParams['font.sans-serif'] = ['SimHei']
# 获取首位的函数
def firstDigital(x):
    x= round(x)
    while x >= 10:
        x //= 10
    return x
# 首位概率累加
def addDigit(lst, digit):
    lst[digit-1]+=1
    return lst
# 理论值：每位概率理论值用于对比
th_freq=[math.log((x+1)/x, 10) for x in range(1,10)]
# 获得2019年报数据
df= ts.get_report_data(2019, 4)
# 只取净利润>0的数据，首先进行次数统计
freq= reduce(addDigit, map(firstDigital, filter(lambda x:x>0, df['net_profits'])), [0]*9)
# 再计算实际概率
pr_freq= [x/sum(freq) for x in freq]
print(th_freq)
print(pr_freq)
# 作图
plt.title('用上市公司2019年报净利润数据验证本福特定律')
plt.xlabel("首位数字")
plt.ylabel("概率")
plt.xticks(range(9), range(1,10))
plt.plot(pr_freq,"r-",linewidth=2, label= '实际值')
plt.plot(pr_freq, "go", markersize=5)
plt.plot(th_freq,"b-",linewidth=1, label= '理论值')
plt.grid(True)
plt.legend()
plt.show()

从图形上看，两者拟合度还是比较高的。
据说有些上市公司数据造假就是被用本福特定律查出来的。所以不认真学习的话，造假都造不好。