• Машинное обучение
• Python
• Data Mining

Однажды томным вечером, сидя напротив мелькающей ленты tjournal и попивая ромашковый чай, внезапно обнаружил себя за чтением статьи про советскую лампочку, которая освещала чей-то подъезд уже 80 лет. Да, весьма интересно, но все же я предпочитаю статьи про политику достижения ИИ в игре дум, приключения ракет SpaceX и, в конце концов, — с наибольшим кол-вом просмотров. А какие вообще статьи набирают внушительные рейтинги? Посты размером с твит про какую-то политическую акцию или же талмуды с детальным анализом российской киноиндустрии? Ну что же, тогда самое время расчехлять свой Jupyter notebook и выводить формулу идеальной статьи.

Сбор базы

Для начала работы стоит раздобыть данные. Поскольку TJ не имеет своего API, то и собирать статьи будем старым добрым скрапингом голого html с помощью requests и вытаскиванием из него полезных данных с помощью BeautifulSoup. На выходе мы получаем 39116 статей с 2014 года. Однако стоит отметить, что сама статья записывается в виде html блока со всеми присущими ей кучей тэгов, а не как простой текст. Сам скрипт тут.

Feature Engineering

После часа работы скрипта у нас на руках есть почти 39000 статей, и пришло время вытаскивать из них фичи всех степеней важности. Но сначала импортируем все необходимое:

import json
import datetime
import numpy as np
import pandas as pd
import snowballstemmer
from bs4 import BeautifulSoup
import itertools

from scipy.sparse import csr_matrix, hstack
from sklearn.feature_extraction.text import TfidfTransformer, CountVectorizer
from sklearn.cross_validation import train_test_split
from sklearn.decomposition import LatentDirichletAllocation
from gensim.models import Word2Vec
import lightgbm as lgb
from sklearn.metrics import r2_score, mean_absolute_error

import seaborn as sns
import  matplotlib.pyplot as plt

%matplotlib inline

Читаем наш датасет и смотрим его содержимое:

dataframe = pd.read_json("data/tj_dataset.json")
dataframe.head()

Поскольку из доступных данных для тренировки мы имеем только контент поста и дату публикации, то и большая часть инженеринга фич будет крутиться вокруг текстовых данных. Из самого поста можно выудить два типа фич: структурные и смысловые (да-да, названия придумал сам). К структурным относятся те, которые, не поверите, определяют структуру поста: кол-во ссылок на ютуб, на твиттер, кол-во слов в названии поста, среднее кол-во слов в текстовых блоках и т.п. Все значения, которые не удалось выловить (обычно из-за малого кол-ва текстовых данных), просто заменям на -9999, градиентный бустинг нас поймет.

Начнем с временных фич. И тут все как нельзя просто — вытаскиваем год, месяц, день, часы и минуты, из инженеринга фич тут только день недели:

def get_time(x):
    time = x.split()[-1]
    h, m = time.split(":")
    return int(h)*60 + int(m)

m_keys = {
    'января': 1,
    'февраля': 2,
    'марта': 3,
    'апреля': 4,
    'мая': 5,
    'июня': 6,
    'июля': 7,
    'августа': 8,
    'сентября': 9,
    'октября': 10,
    'ноября': 11,
    'декабря': 12,
}

dataframe['month'] = dataframe.date.apply(lambda x: m_keys[x.split()[1]])
dataframe['time'] = dataframe.date.apply(get_time)
dataframe['day'] = dataframe.date.apply(lambda x: int(x.split()[0]))
dataframe['year'] = dataframe.date.apply(lambda x: int(x.split()[2][:-1]))
dataframe["weekday"] = dataframe[["year", "month", "day"]].apply(lambda x: datetime.date(*x).weekday(), axis=1)

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

Для начала стоит вытащить из html блока голые куски текста, токенизировать их и пропустить через стеммер полученные слова:

stemmer = snowballstemmer.RussianStemmer()
rus_chars = set([chr(i) for i in range(1072, 1104)])

def tokenize_sent(sent):
    sent = sent.lower()
    sent = "".join([i if i in rus_chars else " " for i in sent])
    words = stemmer.stemWords(sent.split())
    return words

Недолго думая, к списку фич добавляем кол-во слов в посте и заголовке, кол-во html тэгов, среднее и среднеквадратичное отколнение кол-ва слов в блоках. Посты на tj пестрят картинками, видео, референсами на твиттер и ссылками на другие ресурсы, что также можно использовать:

def get_p_data(article):
    html = BeautifulSoup(article, "lxml")
    data = html.findAll("p")
    return [i.text for i in data]

def get_p_count(article):
    html = BeautifulSoup(article, "lxml")
    data = html.findAll("p")
    return len(data)

def get_p_sizes(article):
    html = BeautifulSoup(article, "lxml")
    data = html.findAll("p")
    return [len(preprocess_sent(i.text)) for i in data]

def get_tags_count(article):
    html = BeautifulSoup(article, "lxml")
    data = html.findAll()
    return len(data)

def mean_p(x):
    if not len(x):
        return -9999
    return np.mean(x)

def std_p(x):
    if not len(x):
        return -9999
    return np.std(x)

dataframe['p_list'] = dataframe.article.apply(get_p_data)
dataframe['text_chained'] = dataframe.p_list.apply(lambda x: "\n".join(x))

dataframe['p_list_tokenized'] = dataframe.p_list.apply(lambda x: [tokenize_sent(i) for i in x])
dataframe['p_list_tokenized_joined'] = dataframe.p_list_tokenized.apply(lambda x: list(itertools.chain(*x)))
dataframe['title_tokenized'] = dataframe.title.apply(tokenize_sent)

dataframe['images_count'] = dataframe.article.apply(lambda x: x.count("wrapper-image"))
dataframe['wide_labels_count'] = dataframe.article.apply(lambda x: x.count("wrapper-wide"))
dataframe['link_widget_count'] = dataframe.article.apply(lambda x: x.count("link-widget"))
dataframe['links_count'] = dataframe.article.apply(lambda x: x.count("a href")) - \
    (dataframe.images_count + dataframe.link_widget_count)
dataframe['youtube_count'] = dataframe.article.apply(lambda x: x.count("wrapper-video"))
dataframe['tweets_count'] = dataframe.article.apply(lambda x: x.count("wrapper-tweet"))
dataframe["tags_count"] = dataframe.article.apply(get_tags_count)
dataframe["p_count"] = dataframe.p_list.apply(len)
dataframe["text_sizes"] = dataframe.p_list_tokenized.apply(lambda x: [len(i) for i in x])
dataframe["text_sizes_mean"] = dataframe.text_sizes.apply(mean_p)
dataframe["text_sizes_std"] = dataframe.text_sizes.apply(std_p)
dataframe["text_words_count"] = dataframe.p_list_tokenized.apply(lambda x: len(list(itertools.chain(*x))))
dataframe["title_words_count"] = dataframe.title_tokenized.apply(lambda x: len(x))

Вот с семантическими фичами дело обстоит сложнее. И первое, что приходит на ум: tf-idf — классика в обработке текстовых данных. Но один лишь tf-idf не внушает доверия, так что подключим еще Latent Dirihlet Allocation (LDA) для уменьшения размерности поста до 10-50-и элементного вектора. На всякий случай попробуем взять три LDA по 10, 20 и 50 топиков соответственно.

count_vec = CountVectorizer(tokenizer = tokenize_sent, min_df=10, max_df=0.95)
text_count_vec = count_vec.fit_transform(dataframe.text_chained)
text_tfidf_vec = TfidfTransformer().fit_transform(text_count_vec).toarray()

lda_features = []
topic_counts = [10, 20, 50]
for topics in topic_counts:
    lda = LatentDirichletAllocation(topics, n_jobs=7, learning_method="batch")
    feats = lda.fit_transform(text_count_vec)
    lda_features.append(feats)

lda_features = np.concatenate(lda_features, axis=1)

Можно ли что-то еще интересное добавить? Разумеется, word2vec:

w2v = Word2Vec(dataframe.p_list_tokenized_joined, size=300, workers=7, min_count=5)
w2v_dict = {key:w2v.wv.syn0[val.index] for key, val in w2v.wv.vocab.items()}

vectors = list(w2v_dict.values())
words = list(w2v_dict.keys())
words_set = set(words)

На выходе у нас есть куча векторов слов по 300 элементов в каждом. Теперь нужно их как-то использовать для векторизации поста. Для этого можно использовать 2 техники: вычисление среднего вектора среди всех векторов слов и более сложная, с использованием кластеров.

С первой все понятно — есть матрица c векторами слов, просто берем среднее значение по первой оси, алгоритм обработки пустых постов примерно такой же, как выше :

def get_mean_vec(sent):
    matrix = [w2v_dict[i] for i in sent if i in words_set]
    if not matrix:
        return np.full((300,), -9999)
    return np.mean(matrix, axis=0)

Со второй нужно построить некоторое кол-во кластеров по всем векторам слов и затем присвоить каждому слову в тексте принадлежность к тому или иному кластеру, т.е. получится что-то вроде CountVectorizer-а на стероидах. Так же, как и в случае с LDA, попробуем взять несколько KMeans с разным кол-вом кластеров:

kmeans = []
clusters_counts = [10, 50, 100, 500, 1000]
for clusters in clusters_counts:
    kmeans.append(KMeans(clusters, precompute_distances=True, n_jobs=7).fit(vectors))

word_mappings = []
for kmean in kmeans:
    labels = kmean.labels_
    word_mappings.append(dict(zip(words, labels)))

def get_centroids_vec(sent):
    words = [i for i in sent if i in words_set]
    if not words:
        return np.full((sum(clusters_counts),), -9999)
    result_total = np.asarray([])
    for cnt, mapper in zip(clusters_counts, word_mappings):
        result = np.zeros((cnt, ))
        for word in words:
            result[mapper[word]] += 1
        result_total = np.concatenate([result_total, result])
    return result_total

mean_texts = np.asarray([get_mean_vec(s) for s in dataframe.p_list_tokenized_joined.tolist()])
mean_titles = np.asarray([get_mean_vec(s) for s in dataframe.title_tokenized.tolist()])
clusters_texts = np.asarray([get_centroids_vec(s) for s in dataframe.p_list_tokenized_joined.tolist()])
clusters_titles = np.asarray([get_centroids_vec(s) for s in dataframe.title_tokenized.tolist()])

Все фичи собрали, осталось только скомпоновать их вместе в одну большую матрицу:

feats = ["time", "day", "weekday", "month",
         "images_count", "wide_labels_count",
         "link_widget_count", "links_count",
         "youtube_count", "tweets_count",
         "tags_count", "p_count", "text_sizes_mean",
         "text_sizes_std", "text_words_count",
         "title_words_count"]
X = dataframe[feats].as_matrix()

preprocessed_df = np.concatenate([X, mean_texts, mean_titles, lda_features, clusters_texts, clusters_titles], axis=1)

full_df = hstack([csr_matrix(preprocessed_df), text_tfidf_vec]).tocsr()

Тренировка модели

Поскольку у нас здесь есть кол-во просмотров и рейтинг, то тут прям напрашивается задача регрессии. С каждым месяцем количество подписчиков на tj возрастает, а это значит, что возрастает и кол-во просмотров. Т.е. нужно как-то нормализовать кол-во просмотров, при этом минимизируя влияние выбросов. Для нормализации одного поста можно взять среднее и среднеквадратичное отклонение из предыдущих 200 постов, скалировать текущее значение, при этом убрав все, что находится за первой и 99-ой перцентилями:

hits_y = dataframe.hits

step = 200
new_hits = []
for i in range(step, dataframe.shape[0]):
    hits = hits_y.iloc[i-step: i+1]
    perc_high = np.percentile(hits, 99)
    perc_low = np.percentile(hits, 1)
    hits = hits[(hits > perc_low)&(hits < perc_high)]
    mean = hits.mean()
    std = hits.std()
    val = max(min(hits_y.iloc[i], perc_high), perc_low)
    new_hits.append((val - mean)/std)

hits_y = np.asarray(new_hits)
full_df = full_df[step:]

Глянем на распределение наших рейтингов:

plt.hist(hits_y, bins=100)

Вполне ожидаемый скос налево. Лечится такое обычно логарифмом:

plt.hist(np.log(hits_y - hits_y.min()+1e-8), bins=100)

Уже намного лучше.

min_val = hits_y.min()
hits_log = np.log(hits_y - min_val+1e-6)

Самое сложное позади, осталось только выбрать модель, поставить ее на обучение и молиться на правильный выбор гиперпараметров. Обычно над выбором модели долго не заморачиваются и сразу берут xgboost, что имеет смысл, ибо он уже завоевал свою популярность на kaggle и повсеместно используется как новичками, так и гуру. Но однажды, случайно наткнувшись на сравнение бенчмарков двух реализаций градиентного бустинга, меня привлек внимание LigthGBM от Microsoft своими аппетитными результатами тестов, которые и заставили отложить старичка xgboost как запасной вариант.

В целом, тренировка проходит так же, как и на xgboost — тыкаем туда-сюда max_depth и регуляризацию для избежания оверфита, подбираем наилучшее кол-во деревьев под допустимо минимальный learning_rate.

Датасет для LigthGBM необходимо для начала привести в местный формат:

features = (feats +
    list(itertools.chain(*[[f"lda_{lda}_{i}" for i in range(lda)] for lda in topic_counts])) + 
    [f"mean_matrix_body_{i}" for i in range(300)] + 
    [f"mean_matrix_title_{i}" for i in range(300)] + 
    list(itertools.chain(*[[f"clusters_body_{i}_{j}" for j in range(i)] for i in clusters_counts])) +
    list(itertools.chain(*[[f"clusters_title_{i}_{j}" for j in range(i)] for i in clusters_counts])) +
    [f"tf_idf_{i}" for i in range(text_tfidf_vec.shape[1])])

X_train, X_test, y_train, y_test = train_test_split(full_df, hits_log, test_size=0.2, random_state=753)

train_data = lgb.Dataset(X_train, label=y_train, feature_name=features)
test_data = train_data.create_valid(X_test, label=y_test)

Параметры градиентного бустинга (получены после долгих часов тренировок):

param = {'num_trees': 100000,
         'application':'regression',
         'learning_rate': 0.01, 'num_threads': 7, 
         'max_depth': 10,
         'lambda_l2': 1e-3}
param['metric'] = 'mae'

Непосредственно сама тренировка:

bst = lgb.train(param, train_data, param['num_trees'], valid_sets=[test_data], early_stopping_rounds=200)

Важной особенностью деревьев является то, что есть возможность оценить важность фич:

importance = sorted(zip(features, bst.feature_importance()), key=lambda x: x[1], reverse=True)
for imp in importance[:20]: print("Feature '{}', importance={}".format(*imp))

Feature 'tags_count', importance=546
Feature 'month', importance=510
Feature 'weekday', importance=431
Feature 'images_count', importance=377
Feature 'lda_20_6', importance=332
Feature 'time', importance=316
Feature 'wide_labels_count', importance=291
Feature 'lda_50_48', importance=227
Feature 'text_sizes_std', importance=215
Feature 'mean_matrix_title_134', importance=174
Feature 'title_words_count', importance=161
Feature 'text_sizes_mean', importance=160
Feature 'lda_20_15', importance=149
Feature 'mean_matrix_body_262', importance=142
Feature 'mean_matrix_body_194', importance=141
Feature 'mean_matrix_title_115', importance=138
Feature 'mean_matrix_body_61', importance=136
Feature 'mean_matrix_title_13', importance=136
Feature 'tf_idf_1427', importance=132
Feature 'lda_10_1', importance=130

importance = pd.DataFrame([{"imp": imp, "feat": feat} for feat, imp in importance])
sns.barplot(x="imp", y="feat", data=importance.iloc[:20])

И вот наши результаты:

predictions = bst.predict(X_test)
print("r2 score = {}".format(r2_score(y_test, predictions)))
print("mae error = {}".format(mean_absolute_error(y_test, predictions)))

r2 score = 0.3166155559972065
mae error = 0.4209721870443455

Итоги

C такими параметрами и фичами точность по r2 score получилась 0.317, а по absolute mean error — 0.42. Хороший ли это результат? Вполне неплохо даже с учетом того, что существуют еще пути для улучшения модели. Например, попробовать нейросети (как реккурентные, так и сверточные, с Embedding слоем и весами word2vec), к LDA добавить еще PCA, NMF и прочую декомпозицию, в конце концов, попробовать выделить другие фичи из разряда "отношение кол-ва пунктуации к кол-ву переносов строк". Но дальнейшие изыскания уже оставляю вам.

UPD.
Тут пользователь ben_yazi указал в комментах ссылку на описание апи для TJ.