2013年02月

ユーザー継続率をフィッティング

Rを使って、ソーシャルゲームのユーザー継続率を、カーブフィッティングしてみましょう。 カーブフィッティングとは、データにもっともよく当てはまる数式を見つけ、データに当てはまる曲線を描くことです。 ここで継続率とは、ある日インストールしたユーザーが、n日後にゲームを訪問している確率と考えます。 たとえば、以下はある日にインストールしたユーザーが日を経過するにつれてへっていくところを10日分集計したデータです(ただし数字はダミーです)。
installがインストール日に加入したユーザー数、nが経過日数、countはインストールユーザーの内、n日後にも継続してアクセスしてきたユーザー数です。最初は1301人いたユーザーが10日後には、168まで減っています。


installncount
13011338
13012264
13013248
13014235
13015185
13016173
13017181
13018188
13019145
130110168



上記のようなデータをインストール日5日分集め、データからえられた継続率(count/install)をグラフにプロットしたのが以下の図です。
retention

ここでは、ユーザーの期待値を求めたり、未来のユーザー数を予測するために、この減少率を数式で表現し、曲線でフィッティングしてみたいと思います。
先に結果を出しておくと、以下が継続率にフィッティングした曲線です。
retention-fitting

様々なインストール日のデータに当てはまるような曲線が引けていると思います。当然ながらここから定着するユーザー数を予測することもできます。
以下では、この曲線がどのような式に基くものか解説していきましょう。


考え方

継続率のような確率的なデータを扱うには、通常の回帰分析ではなく、一般化線形モデルの一種であるロジスティック回帰モデルという統計モデルを利用します。
このモデルに基づくと、n日後の継続率は以下のような式によって表わされます。



logisticという関数は以下のような形の関数で、値域が(0,1)となる関数です。確率が0未満や1より大きくなってしまっては困るので、この関数を利用して(0,1)区間にマッピングしていると考えてください。



また日数nのlogをとっていますが、これは継続率がしばしば曲線的な形状を取るため、こちらの方が当てはまりがよくなるという経験則に従っているだけです。サービスの性質によっては必要ないかもしれません。

問題は、(1)式のαとβを調整することによって、過去のデータにもっともよくはてまるα、βの値を見つけることです。より正確に言うと、インストールしたユーザー数に対する継続したユーザー数を確率P(n)に従う二項分布と見なした上で、もっともよく当てはまるα、βを見つけます。一般的に用いられるのは最尤法と呼ばれる手法です。最尤法では、尤度、つまり与えられたデータが得られる確率を最も高くするα、βの値を探します(この辺の解説は最後に書籍を紹介するのでそちらを参考にしてください)。

使用するデータは以下のようになっています。上にあげたデータと同様、ある日にインストールしたユーザーを10日分追跡したデータを、インストール日5日分集めたものです。

dummy <- read.table("dummy.csv",sep=",", head=T)

> dummy[1:5,]

   install  n count

1     1301  1   338

2     1301  2   264

3     1301  3   248

4     1301  4   235

5     1301  5   185


一般化線形モデルはRがデフォルトでサポートしているのでそちらを利用します。glm関数のfamily引数にbinomial(二項分布)を指定します。

inputdata <- cbind(dummy$count, dummy$install-dummy$count)
n <- dummy$n
model <- glm(inputdata~log(n), family=binomial)

最初の二行は、データの準備を行なっているだけで、メインはglm関数を使用している一行だけです。
~(チルダ)の左側が被説明項で、ここでは残存ユーザー数と離脱ユーザー数を水平に結合した行列になっています。~の右側は説明項です。familyにbinomialを指定すれば、先程説明した最尤法によるパラメーターの探索を行ない、モデルを構築してくれます。
モデルを構築するにはこれだけで十分です。
係数を表示するには、model$coeffを見てください。

> model$coeff

(Intercept)      log(n)

-0.9853399  -0.4220020


(Intercept)とあるのがグラフの切片α、log(n)がlog(n)の係数βです。
つまり、このデータからえられるユーザーの継続率は以下の数式で表現できます。
これによって、過去の継続率に基づいて、任意の日数後の継続率を求めることができました。



Rで継続率の予測値をえるには、以下のようにします。

logistic <- function(x) return(1/(1+exp(-x)))

coeff <- model$coeff

n <- seq(1,15)

# 継続率の予測値

yy <- logistic(coeff[1] + coeff[2] * log(n))


# グラフにプロットしてみる

plot(dummy$n,dummy$count/dummy$install,xlab="day",ylab="retention rate")

lines(n, yy, col=3)

retention-fitting


書籍紹介
 
一般化線形モデルについては以下の本が大変わかりやすく説明も豊富ですのでおすすめです。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

 

pythonでElastic MapReduce

Geographica restituta per globi trientes

まとめ

Elastic MapReduceは、amazonのAWS上でHadoopを使えるサービスです。一時的にインスタンスをたくさん立ち上げることで、重めのバッチ処理を分散処理できます。料金はEC2インスタンス分 +αで使用できます。
楽! インスタンスいっぱいたちあげます、計算します、結果まとめますという一連の処理を気軽に書けて、コマンド一発で実行できます。特に集計処理、バッチ処理には便利です。
ただデバッグは手間がかかります。


MapReduceについて

MapReduceは、大量のデータを複数のマシンで分散して扱うための技術です(デザインパターン的な)。
基本的な考え方は、処理を以下の2種類に分けるだけです。
  •  map: データ1行に対する処理
  •  reduce: 集計処理



mrjobの紹介

Hadoopは基本的にJavaで書くんですが、Hadoop Streamingと言って他の言語の標準入出力を利用することもできます。今回はこちらを利用してます。
pythonのmrjobというライブラリをつかえば、Hadoopのめんどうな部分をラップして使うことができます。

動かすための設定ファイルは以下のようになります。個々のパラメーターは実行時オプションでも上書きできますが、重要なのはAWSにアクセスするための aws_access_key_id と aws_secret_access_key です(アクセスキー自体はAWSのWebコンソールで発行できます)。
↓以下のファイルを ~/.mrjob.confという名前で保存
 
runners:
  emr:
    ec2_instance_type: m1.large
    num_ec2_instances: 4
    aws_access_key_id: ひみつ
    aws_secret_access_key: ひみつ


実行するスクリプト

サンプルとして、アクセスログから時間帯ごとのPVを抽出するスクリプトを掲載します。
#!/usr/bin/env python
# coding:utf-8

import traceback
from cStringIO import StringIO
import csv
import datetime
import mrjob
from mrjob.job import MRJob
import sys

class CountHourly(MRJob):
    def configure_options(self):
        super(CountHourly, self).configure_options()
    def __init__(self, args=None):
        super(CountHourly, self).__init__(args)
    def mapper(self, key, value):
        # mapper: ログの各行を処理する部分。ここではアクセスログをパースし、時間を取得している。
        try:
            reader = csv.reader(StringIO(value), delimiter=' ', quotechar='"')
            row = reader.next()
            if not row:
                return
            actime = datetime.datetime.strptime(row[1]+' '+row[2], '%Y-%m-%d %H:%M:%S')
            hour   = actime.hour
            key = [hour]
            yield hour, 1
        except Exception:
            sys.stderr.write(traceback.format_exc() + "\n")

    def reducer(self, key, data):
        # reducer: mapperのアウトプットを集計する。ここでは時間帯別に データの長さを取得するだけ。
        yield key, sum(data)

if __name__ == '__main__':
    CountHourly.run()



以下のように実行します。はじめはローカルでデバッグするのがおすすめです。


# デバッグ用にローカルで
python counthourly.py path/to/log > result
# elastic mapreduce利用
python counthourly.py -r emr s3://mybucket/path/to/log > result


ディレクトリを引数に指定すると、その下の全ファイルを入力として実行します。またbz2形式は勝手に解凍して処理してくれます。
内部でインスタンスを起動し、スクリプトをコピーし、実行して結果をS3から取得しますが、そのあたりはmrjobがラップしてくれます。

データ分析グループのブログを立ち上げました

Data Center
はじめまして、KLab株式会社のデータ分析グループです。
ソーシャルゲームのデータ抽出やデータマイニングのシステムをつくっている開発者のチームです。 
ソーシャルゲームのデータマイニングにかかわる技術やツールの紹介をしていこうと思います。