Закон Бэнфорда и COVID-19

Утро вторника, чашка кофе, интенсивная борьба с выгоранием, бумажная работа, на глаза попадается датасет с данными по COVID-19. «Ого, это интересная штука, много цифр по какой-то тематике, что с этим обычно делают?», подумал я.

Поисковик предложил эти данные проанализировать, раз я не программист, то анализировать питоном, раз я не разбираюсь в теории вероятности, то законом Бэнфорда. Любопытство: «Ок, бро, ты в деле!».

Закон Бэнфорда подошел мне по нескольким причинам. Датасет имеет числа не ограниченные диапазонами значений, достаточно большой, содержит цифры и эпидемии имеют экспоненциальный рост.

Python подошел мне по нескольким причинам. Я люблю змей. Два года пытался его изучать в теории, но оказалось это совершенно не нужно, достаточно просто начать программировать на нем.

Датасет представлял собой три варианта хранения, я выбрал json, потому что мне нравится как это звучит и я слышал, что это крутая и опасная штука в умелых руках.

Для начала я сдернул 24 вкуснейших мегабайта к себе на диск, чтобы на время тестов не дергать каждый раз это с сервера:

import json
import urllib.request


def get_json(url_data):
    web_url = urllib.request.urlopen(url_data)
    data = web_url.read()
    data = json.loads(data)
    return data


def corona_get():
    url = 'https://covid.ourworldindata.org/data/owid-covid-data.json'
    json_data = get_json(url)
    with open('data.json', 'w') as f:
        json.dump(json_data, f)

Затем, как фильме Нолана «Начало» я спускался все глубже в недра json, пока не достиг нужных данных total_cases. Как оказалось, датасет имеет не очень обязательную структуру.

Во избежание исключений я сделал проверку через try-except собирая число ошибок, для определения погрешности (на самом деле было бы неплохо где-то выводить число ошибок на срез, но меня захватила сама идея и я просто заленился).

Затем перешел непосредственно к алгоритму определения частоты появления цифры.

        
# Роемся в данных на регион, при желании можно проверять документ не вдоль, а поперек, делая срезы по датам, но об этом позже.
for e in json_data[region]['data']:
    # Так я удумал проверять наличие значения в документе json.
    def data_cases_exist():
        try:
            test = e['total_cases']
            return True
        except:
            return False
    # Если значение существует, таки добавляет его в список. 
    if data_cases_exist():
        list_cases.append(int(e['total_cases']))
    # Иначе считает сколько таких ошибок на срез, неплохо бы такое знать, вообще-то.
    else:
        error_count += 1

Уберите ребенка от экрана, потому что далее идет код, требующий еще большего рефакторинга, так как частично является заимствованием с неуклюжей доработкой.

В нем я создаю отсортированный для распределения по закону Бэнфорда словарь, где ключом является первый разряд числа, а значением частота его появления.

# умные люди говорят, что с коллекциями проще, чем со словарями.
import collections    

# Делает лист benford_list из первых цифр значений
for i in list(map(str, list_cases)):
    benford_list.append(i[0])

# Тут питон делает подсчет частоты появления цифры в benford_list
results = collections.Counter()
for i in benford_list:
    results1 = collections.Counter(str(i))
    results = results + results1

# Сортировка коллекции в порядке убывания по ключам, словарь якобы так красиво не умеет. "Подписывайесь на канал, ставьте лайки и оставляйте коментарий так ли это."
dict_b = dict(results.most_common(10))

# Тут я вроде бы пытался конвертировать строковые ключи в целое, зачем, уже не помню сам. 
dict_b = {int(k): int(v) for k, v in dict_b.items()}
    
# Сортровка по ключам от 1 до 9
sorted_dict = dict(sorted(dict_b.items()))

# Ноли нам не понадобятся, поэтому убираем ноль из нашего словаря
sorted_dict.pop(0)
return sorted_dict

И вот, наслаждаясь своей работой, я держу в руках голубоватый, переливающийся на свету отсортированный sorted_dict региона RUS.

{1: 14559, 2: 10128, 3: 3798, 4: 4431, 5: 3587, 6: 4009, 7: 3798, 8: 4853, 9: 4853}

Показываю друзьям. Кидаюсь ссылкой на numberphile о законе Бэнфорда. Рассказываю про эпидемии… Что ж, хвастаться числовыми словарями перед друзьями получилось не очень, придется визуализировать красивыми гистограммами.

Окей гугл, чем строят графики хорошие мальчики в 2020? Гугл голосом {username} отвечает: «Есть одна очень крутая штука, может всё, если хорошо разберешься». Если хорошо разберешься… Challenge accepted!

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

# Передаем наш словарь и наименование региона, чтобы сразу понимать, что на гистограмме.
def benford_plot(sorted_dict, region):

    # x это наши цифры от 1 до 9 
    x = []

    # y1 это частота появления цифры
    y1 = []

    # y2 это список по которому мы будем строить кривую по закону Бэнфорда для сравнения.
    # Насколько я понимаю трактовку закона, танцевать нужно от 1, либо... нуакак?
    y2 = [sorted_dict[1]]
    
    # Строим остальные 8 столбцов по закону Бэнфорда, первый, по идее, это максимум.
    for i in range(8):

        # Каждый следующий стобец ниже на ben_percent
        ben_percent = math.log10(1 + (1 / (i+1)))
        value = y2[i] - (y2[i] * ben_percent)
        y2.append(value)

    # Наполняем листы для отображения осей координат 
    for key, value in sorted_dict.items():
        x.append(key)
        y1.append(value)

    fig, ax = plt.subplots()
    
    # Гистограмма по данным из датасета
    ax.bar(x, y1, color=[0, 0, 0], width=0.5)
    ax.xaxis.set_major_locator(ticker.MultipleLocator(1))

    # Гистограмма по мнению товарища Бэнфорда
    ax.bar(x, y2, color=[0.4, 0.4, 0.4, 0.35415316])

    # Расставляет красиво значения над столбиками
    for x, y in zip(x, y1):
        plt.text(x, y + 0.05, '%d' % y, ha='center', va='bottom')

    fig.set_figwidth(12)
    fig.set_figheight(6)
    fig.set_facecolor('floralwhite')
    ax.set_facecolor('seashell')

    plt.title('Total cases in ' + region + ' for all time', fontdict=None, loc='center')
    plt.savefig('pic\\' + region + '.png')

    plt.close()

Таким образом мы сложили в pic гистограмму с регионом и если определить все это в цикл, перебирающий в json все регионы, то получим еще больше красивых картинок. Поэтому:

if __name__ == '__main__':
    f = open(r'data.json')
    json_data = json.load(f)
    for i in json_data:
        location = json_data[i]['location']
        benford_plot(benford_data(i), location)

Для проверки я использовал генератор случайных чисел и ряд Фибоначчи, потому что первый есть хаос, а второй есть экспонента, ровно как и рост заболеваемости при эпидемиях.

Иными словами ГСЧ не впишется в закон Бэнфорда, а ряд Фибоначчи и все графики с датасета впишутся идеально.

Для контроля статистической значимости я использовал Z-тест, ведя расчет по формуле: $inline$\frac {|AP - EP|-\frac{1}{2N}} { \sqrt {\frac { EP(1-EP)} {N} }}$inline$, взятой из нескольких источников, например отсюда или отсюда, где AP вероятность распределения из датасета, EP теоретическая вероятность по Бэнфорду, N число иcследуемых элементов. Исходя из доверительной вероятности в 1,96 на графиках внизу столбца Z < 1,96 выводил белым, Z > 1,96 оранжевым и Z > 10 красным. Запутались? Я запутался, если вы понимаете какие выводы можно сделать из этих значений, пишите, попробуем разобраться. Насколько я понимаю, при достаточно больших данных такое возможно, что z-stest показывает все значения как статистически значимые и принимать во внимание нужно действительно высокие показатели Z, именно поэтому я выбрал градацию цветов с оранжевым и красным. Говоря проще, столбцы с оранжевыми показателями настораживают, а с красными вызывают явные подозрения в нарушении закона Бэнфорда.

Закон Бэнфорда и ряд Фибоначчи (50000 чисел):



Закон Бэнфорда и случайные числа (50000 чисел):




Такие вот интересные результаты. За время написания материала я нашел ошибку в алгоритме, переписал код, нашел ошибку в коде, опять переписал, а еще меня посетило много мыслей, которые по-хорошему должны посещать до написания статьи. И я подумал, а может я где-то ошибся? Для чистоты эксперимента я попробую найти датасет по эпидемиям, данные о которых с меньшей вероятностью кто-то мог пытаться фальсифицировать, например туберкулёз, но на поверку такие данные найти не так и просто.


Вдохновлялся Testing Benford’s law.

Данные c ourworldindata.org.