読者です 読者をやめる 読者になる 読者になる

広告/統計/アニメ/映画 等に関するブログ

広告/統計/アニメ/映画 等に関するブログ

バナー広告の接触回数分布を関数で回帰してみる

広告業界ではある媒体に出稿した時の接触回数の分布(1回当たった人がn人、2回当たった人がn人)をベータ二項分布などで近似しシミュレーションするのが一般的と言われている

roi-plus.com

なおこのサイトによるとWeb広告の場合は、負の二項分布がグローバルで標準らしい。*1

ということで、バナー広告の接触回数分布がどんな関数になるのか、入手可能なデータで確かめてみた。

ビデオリサーチ社が2年前に調査したデータを使う

データソースについて

説明はここを参照して欲しい。

hajimeteweb.jp

リリースの最後の方に、接触回数の分布がある。

プレスリリース 2015年 | ビデオリサーチ

このサンプル数がわからなかったが、計算の簡単のためこのパーセンテージの小数点が消えるように桁を変えて,仮にこの値を回帰できる関数を探してみる。

x    y
1   271
2   179
3   129
4   89
5   69
6   53
7   40
8   32
9   25
10  20

プロットしてみる

プロットした時の関数は省略する。

一般的に言われている負の二項分布っぽい感じはする。

f:id:yyhhyy:20170416184626p:plain

負の二項分布について

負の二項分布の理解

負の二項分布は、確率pの試行をN回繰り返す時に、r回成功するまでに何回試行を繰り返すか?の分布を知る為に使われるのだが、今回の場合はその理解というよりも、こちらのブログの通り、

sugisugirrr.hatenablog.com

カウントデータの回帰にポワソン分布を使おうとしたら分散がもっと大きい分散なので、ガンマ分布とポアソン分布の階層モデルを作るとき、それが負の二項分布になる、というポアソン分布の過分散時の推定、と言われている方の用途で使う。

負の二項分布について

このサイトを見ると確かにガンマ分布で何パターンか作ったポアソン分布を合成していくと負の二項分布になっている。。。

d.hatena.ne.jp

負の二項分布 | r回の成功を得るのに必要な試行回数

RStanに負の二項分布のあてはめ

Rのパッケージ、nbGLMでやってもほぼおなじ結果なのだがせっかくなのでRStanでベイズモデリングしてみる。

知識不足のため、収束させる為の条件はわからなかったのでコードはほぼほぼこれらのサイトを見よう見まねで使った。

Simulate and fit negative binomial GLMs in Stan

Stanコードの書き方 中級編

モデリングのコード

接触回数xを入れたら、接触人数yが予測できるようなモデリングを仮定する。線形回帰と異なってリンク関数なので、y直接ではなく、yの自然対数を回帰する。

log(y) = b0 + b1 * x

モデルのコードは別ファイルにしてもいいが、個人的には同じファイル無いで「"“」で挟んで代入する方が見直す時に便利なので好き。

model <- "
data {
  int<lower=0> N;
  vector[N] x;
  int<lower=0> y[N];
}
parameters {
  real<lower=0> phi;
  real b0;
  real b1;
}
model {
  phi ~ cauchy(0,3);
  b0 ~ normal(0,5);
  b1 ~ normal(0,100);
  y ~ neg_binomial_2_log(b0 + b1 * x,phi);
}
"

データの読み込みとRStan用のリスト化

RStanに受け渡す為にデータをリスト化する

d <- read.csv("banner.csv",header=T)
data <- list(N=nrow(d),x=d$x,y=d$y)

実行

実行した後、summary(fit)で先程仮定した「b0」や「b1」の平均値や四分位数が出て来るが、meanを使って値を代入すればよい。

実践 ベイズモデリング -解析技法と認知モデル-

実践 ベイズモデリング -解析技法と認知モデル-

のコードを見て、”なるほど”と思った。流石です。

fit <- stan(model_code=model,data=data,seed=1234)
rs <- rstan::extract(fit)
b0 <- mean(rs$b0)
b1 <- mean(rs$b1)

暫く待つと以下の値を得た。*2

> (b0)
[1] 5.794788
> (b1)
[1] -0.2975039

回帰の関数

これで、x(接触回数)からy(接触人数)を予測する関数ができた。log(y)を回帰しているので、ネイピア数に計数を累乗すればいい。

y = exp( 5.794788 +  -0.2975039 * x)

元のデータより広い値を入れてみる

恐らくビデオリサーチのデータは、10回よりも多いデータもあった筈だ。(パーセンテージを足すと100%にならない)

そこで、新しく d_predという名前のデータフレームを作り、そこに先程の関数で予測される値を格納し、そのデータと元のデータをdplyrのレフトジョインで一つのデータフレームに入れてみる

x_fit <- seq(1:20)
y_fit  <- exp(b0 + b1 * x_fit)

d_pred <- as.data.frame(NULL)
d_pred <- as.data.frame(x_fit)
colnames(d_pred) <- c("x")
d_pred$y_fit <- y_fit

d_pred <- dplyr::left_join(d_pred,d)
> head(d_pred)
  x     y_fit   y
1 1 244.02836 271
2 2 181.23247 179
3 3 134.59587 129
4 4  99.96028  89
5 5  74.23747  69
6 6  55.13392  53

余り精度は高くないようだが、これが今の自分の知識の限界なので勘弁して欲しい。

実測値と予測値を並べてグラフにしてみる

実測値と予測値の棒グラフを並べてみる

p <- ggplot(d_pred_m,aes(x=x,y=value,fill=variable))
p <- p + theme_bw(base_family="Japan1GothicBBB")
p <- p + geom_bar(stat="identity",position = "dodge")
p <- p + geom_text(aes(x=x,y=value,label=round(value)),vjust=-0.5,position = position_dodge(1))
p <- p + scale_fill_brewer("Set2")
p <- p + xlab("接触回数")
p <- p + ylab("人数分布")
plot(p)
ggsave(plot=p,file="20170416-1.png",dpi=300,width=4,height=3,scale=2)

個人的な感覚では、一回しか接触しないバナーが殆どだとは思うけれど、2015年のビデオリサーチの大規模調査を信用するとするとこんな感じである。自分は接触しなくても、ヘビーにWebを使う人ととそうでない人とバラつきがあるということだろう。

f:id:yyhhyy:20170416191802p:plain

媒体数が限られているマス広告と比べて、果してバナー広告の接触回数分布が簡単に回帰できる類いのものなのか?という疑問はあるし、そもそもターゲティングが違えばまた違うだろうし、そもそもリマーケティングを入れたりフリークエンシーキャップをかけたらもっと複雑になるので予測すること自体が難しいのではないか?とは思う。

ただ、「負の二項分布で近似できる」という話については、嘘ではなさそうだ。

*1:文献に全然ヒットしないけどきっとどこかにあるのだろう

*2:省略するがMASSパッケージの標準の回帰の関数でも同じような数値が出た

アンケートや視聴率の誤差を推測するプログラム(3)-プログラムのスタンドアローンな実行ファイル化-

これは3回中3回めの記事です

Pythonがない人の環境でもプログラムを実行できるように実行ファイル化します。色々調べるとPython3ではまだ最新版に対応されていないようで、Python2の環境を手に入れるところから始まります。

先ず、Python2.7の環境を手に入れる

1からインストールするでも良いですし、default以外に2.7用の環境を作るでもいいです

python3とpython2を共存させる

コマンドプロンプトで以下のコマンドで作るだけでよい

conda create -n py2 python=2.7 anaconda

これで新たにpython2の環境が構築される。切替えや現在のデフォルトを知るには、以下のようにする

#デフォルトの環境がどっちになっているか確認する
conda info -e

#Python2の環境に変項する
activate py27

#デフォルトの環境に戻す
deactivate py27

参考にしたサイト

tug-uca.hatenablog.com

自分のPCはプロキシの問題があってエラーが出てしまったのでAnacondaをPython2.7のAnaconda2にまで下げて以下を実行しています。

pyinstallerを手に入れる

スタンドアローン化するパッケージは幾つかあるようだがpyinstallerが一番手軽。

以下のサイトの夫々のページのアドバイスに沿ってpyinsatallerを導入する

参考にしたサイト

PythonでWindowsアプリケーションを作ってみよう – Nobwak's Lair

Windowsにpyinstallerをインストール – Nobwak's Lair

pywin32が必要

先程のリンクの2つめ。Python2.7用をダウンロードする。インストーラー型で配布されているので、普通に実行すれば良い。

pyinstallerを入れるときの注意

Pythonのパッケージが入っているフォルダはここなので、ここにpyinstallerなどのフォルダを自分で作っていれる。

C:\(自分でインストールしたAnacondaの場所)\Lib\site-packages\

pip installが使えないとき

会社のパソコンなどを使っていると、往々にしてプロキシがありpip install できないことが多い。そういう時はパッケージを格納したディレクトリまでコマンドプロンプトで移動し

python setup.py install

とする

途中で足りないパッケージがあるよとerrorで言われたら、ググってパッケージをダウンロードし、同様にsetup.pyをPythonで実行してインストールしてやり直す。

pyinstallerを作ってみる

.pyのファイルを確認する

コマンドプロンプトで格納したファイルのディレクトリまで移動し

python hogehoge.py

として実行してみる。

元々のプログラムが上手く行っているか行ってないか事前に確認しておく。

pyinstallerを使う

大丈夫であればpyinstallerを使う。

参考にしたサイト

fnwiya.hatenablog.com

torina.top

qiita.com

コマンド自体はとても簡単。

pyinstaller hogehoge.py --onefile --clean
  • 単独のファイルにするには「–onefile」のオプションを。
  • やり直しする場合は、「–clean」で前のファイルを必ず消す。
  • コマンドプロンプトを立ち上げないようにするには、「–noconsole」をつける

なお、実行には10分位は待つ。

そして実行ファイルをダブルクリックしてからかなり待つ

とても便利だけど普通の人にはちょっと待ちすぎ?という位長い時間起動にかかるのがたまにきず。ファイルも重い。

参考にしたサイト

www.slideshare.net

d.hatena.ne.jp

Excelでもできる

ところでここまでやってからExcelの標準関数でできることがわかった。

例えば、サンプルサイズを「C2」のセルに、値(視聴率)を「C3」のセルに入力させると、

=ROUND(C3 + NORMSINV(0.025)*SQRT(C3*(100-C3)/C2),2)
=ROUND(C3 + NORMSINV(0.975)*SQRT(C3*(100-C3)/C2),2)

とすればよい

f:id:yyhhyy:20170408232635p:plain

参考にしたサイト

dekiru.net

アンケートや視聴率の誤差を推測するプログラム(2)-プログラムのGUI化-

このエントリーは3回中の2回目です

Pythonが入っているPCであれば先般の記事のプログラムで充分なのですが、コマンドプロンプトの画面に数字を打ち込むのは一般の人には抵抗があると思われます。

余談ですがコマンドプロンプトの背景は白の方が心が落ち着くのでオススメです!

f:id:yyhhyy:20170408215142p:plain

というわけで、GUIベースのプログラムに書きかえてみます。

Tkinterを使う

基本的な使い方はこちらの通り

www.shido.info

qiita.com

root = Tk.Label(None)

root.mainloop()

と記載してこの2行の間にドンドン追加していきます

入力値を取得して表示する関数を作る

サンプルサイズや結果の値を「samplesize_entyr.get()」や「p_ori_enty.get()」などで取得するとして、取得した値を数値化「float()」し、先程計算した一連の計算をし、表示する関数を作ります。

def kakunin():
    n_samples = float(samplesize_entry.get())
    p_ori = float(p_ori_entry.get())
    p = p_ori/100.0
    
    a,b = st.norm.interval(alpha=alpha, loc=0, scale=1)
    
    gosa = np.array([a,b])*np.sqrt(p * (1 - p) / n_samples)
    
    rs = (p + gosa) * 100
    rs_r = [round(rs[0],2), round(rs[1],2)]
    
    anser = Tk.Label(text = u"標本平均に対する95パーセント信頼区間の値は",font = ("14"))
    anser.grid(row = 4,column =0,columnspan =2)
    
    p_hani = Tk.Label (text = u" %(s1)s パーセントから %(s2)s パーセントの間" % {"s1":rs_r[0],"s2":rs_r[1]},font = ("14"))
    p_hani.grid(row = 5, column = 0,columnspan =2)
小数点の四捨五入

xに対して、小数点以下nまでで四捨五入するには以下の数式を使います

round(x,n)

配列全体には使えないので上のコードでは要素一つ一つにroundを使っています

配列の表示

テキスト分の中に、配列の数値を与えるには、一度データを辞書型に変換してキーを与えるのが便利です

Tk.Label (u" %(s1)s パーセントから %(s2)s パーセントの間" % {"s1":rs_r[0],"s2":rs_r[1]})

なおpython2では日本語は必ずUnicodeの指定「u」が必要です

参考にしたサイト

www.lifewithpython.com

Tkのウィジェットを並べる

Tkinterのgridを使って縦横どこに置きたいかを指定します。例えば、サンプルサイズを入力して欲しい、というコメントと入力欄を並べるには、以下のように書きます

samplesize_label = Tk.Label(text = u"サンプルサイズを入力⇒")
samplesize_label.grid(row= 0,column = 0)

samplesize_entry = Tk.Entry()
samplesize_entry.grid(row = 0, column = 1)

f:id:yyhhyy:20170408221314p:plain

参考にしたサイト

www.shido.info

関数を適用させする

クリックで実行するように、ボタンに対して「command」で先程の関数を適用させます。

main_Button = Tk.Button(text = u"入力したらこのボタンをクリックして下さい",command = kakunin)
main_Button.grid()

参考にしたサイト

stackoverflow.com

コード全体

順番にウィジェットを並べて配置していき、最後に表示させます。

# coding: utf-8
import numpy as np
from scipy import stats as st
import Tkinter as Tk

alpha = 0.95

root = Tk.Label(None)

def kakunin():
    n_samples = float(samplesize_entry.get())
    p_ori = float(p_ori_entry.get())
    p = p_ori/100.0
    
    a,b = st.norm.interval(alpha=alpha, loc=0, scale=1)
    
    gosa = np.array([a,b])*np.sqrt(p * (1 - p) / n_samples)
    
    rs = (p + gosa) * 100
    rs_r = [round(rs[0],2), round(rs[1],2)]
    
    anser = Tk.Label(text = u"標本平均に対する95パーセント信頼区間の値は",font = ("14"))
    anser.grid(row = 4,column =0,columnspan =2)
    
    p_hani = Tk.Label (text = u" %(s1)s パーセントから %(s2)s パーセントの間"                          % {"s1":rs_r[0],"s2":rs_r[1]},font = ("14"))
    p_hani.grid(row = 5, column = 0,columnspan =2)
    
    syuryo = Tk.Label(text= u"右上の×印をクリックして終了")
    syuryo.grid(row = 6, column =0, columnspan =2)

title = Tk.Label(text=u"標本平均の95%信頼区間での標本誤差を推定します",font=("14"))
title.grid(row = 0,column = 0, columnspan = 2)

samplesize_label = Tk.Label(text = u"サンプルサイズを入力⇒")
samplesize_label.grid(row=1,column = 0)

samplesize_entry = Tk.Entry()
samplesize_entry.grid(row = 1, column = 1)

p_ori_label = Tk.Label(text = u"パーセンテージを入力(5.5パーセントなら、5.5と入力)⇒")
p_ori_label.grid(row = 2,column =0)

p_ori_entry = Tk.Entry()
p_ori_entry.grid(row = 2,column = 1)

main_Button = Tk.Button(text = u"入力したらこのボタンをクリックして下さい",command = kakunin)
main_Button.grid(row = 3,column =0,columnspan =2)

root.mainloop()

数値を入力してクリックするとこのような画面で表示されます

f:id:yyhhyy:20170408222133p:plain

アンケートや視聴率の誤差を推測するプログラム(1)-標本誤差の計算まで-

アンケート調査や視聴率の数字を見せた時に、 「これはn数が少ないから信憑性が薄いねぇ」 と言われることがあると思います。しかし実際のところ、 ”では、今のn数だとどれくらいの幅で理解すべきなのですか?” という質問に応えてくれる人は少ないでしょう。というのも意外とその場で計算するのは容易ではないからです。

今回は、サンプルサイズから標本誤差を推定するスタンドアローンな実行ファイルを作るまで、3回に分けてブログにします。

標本誤差について

どんな統計の本にもあるような基礎的な話ですが、国勢調査などを除いて、殆どの調査は全数調査はしません。

サンプルを集めて、該当質問の平均値を取り、それを結果の値とします。

更に、以下にランダムにサンプルが集められていたとしても偏りが発生していないという絶対的な保証はありません。そのような時、 少なくとも平均値については、サンプル数が多ければ大数の法則がなりたつ ということを使って、 同じような調査を何度も何度もやったらそのうち95%は平均的に○○~○○の値に収まるであろう という何とも消極的なことしか分からないのです。

標本誤差とは?

どういう計算をするのか?はこのサイトが分かりやすいでしょう。

www.stat.go.jp

ビデオリサーチ社の視聴率の場合

標本誤差 | ビデオリサーチ

サンプルサイズから結果の値の標準誤差を知る為のプログラムを作る

  • 普段統計の計算はRを使っていますが、今回は最終的に実行ファイルにしたいのでPythonにします。
  • 現時点では実行ファイルにするにはPython2の方が良いのでPython2を使います

Pythonを使った標本誤差の計算

正しい手法はこのサイトの通り

Python正規分布の平均値の信頼区間を計算する方法 」

document

しかし先程の参照サイトの通り、母数に対してサンプルサイズがすくない場合はt分布ではなく正規分布でも良い、とありましたので、簡略化することにします。

Pythonのscipy statsを使う

使うモジュールは、scipyの中のstatsから、norm.intervalという関数を使って、平均値0、分散1の正規分布の時の95%の有意点の上下を求めます

参考にしたサイト

kaisk.hatenadiary.com

from scipy import stats as st
alpha = 0.95
a,b = st.norm.interval(alpha=alpha, loc=0, scale=1)

すると、以下のような数値が出てきます。ビデオリサーチ社のサイトでも「2」を使っていましたが、妥当な値と言えるでしょう。

(-1.959963984540054, 1.959963984540054)

サンプルサイズやアンケート結果の値を画面から入力する

サンプルサイズをコード本文にベタ打ちすると、毎度毎度コードの修正が必要になってしまうので、コンソールから都度都度入力するようにしましょう

参考にしたサイト

qiita.com

print(u"サンプル数を入力して下さい")
n_samples = float(input('>>>  '))
print (u"パーセンテージを入力して下さい。(5.5パーセントなら、5.5と入力)")
p_ori = float(input('>>>  '))

このような感じで自分で数字を打込み、それが引数となるようになります

f:id:yyhhyy:20170408213409p:plain

誤差の値を計算する

一番最初の参照サイトの簡略化された数式を実行します。配列のためだけにnumpyを使っていますが、計算自体はとてもシンプルです。

import numpy as np
p = p_ori/100.0
gosa = np.array([a,b])*np.sqrt(p * (1 - p) / n_samples)

仮にビデオリサーチ社のサイトのように900サンプルで視聴率10%だったとして入力すると、サイトの通り、プラスマイナス2%位は割り引いて考える必要があります。*1

[-0.01959964  0.01959964]

コード全体

今回のコード全体は以下のようになります

# -*- coding: utf-8 -*-
import numpy as np
from scipy import stats as st
alpha = 0.95
print(u"サンプル数を入力して下さい")
n_samples = float(input('>>>  '))
print (u"パーセンテージを入力して下さい。(5.5パーセントなら、5.5と入力)")
p_ori = float(input('>>>  '))
p = p_ori/100.0
a,b = st.norm.interval(alpha=alpha, loc=0, scale=1)
gosa = np.array([a,b])*np.sqrt(p * (1 - p) / n_samples)
print (u"標本平均の95パーセント信頼区間での値は次の間")
print ( (p + gosa) * 100)
print (u"単位パーセント")

*1:とてもややこしい話ですが、95%の確率で視聴率は8-12%だ!というわけではないのが伝統的な統計学の難しさです。同じような調査を100回やったうち95回は、8-12%の間であろう、というだけで真の値がどこかはわからないのです。

【備忘録】Rでディレクトリ配下の複数ファイルからデータをデータフレームに格納する

ディレクトリのファイルパスを指定する

例えば、カレントディレクトリの下に「data」というフォルダを作ってその中から読ませたければ、

path <- 'data/'

などとしておく

※参考にしたサイト

R:list.filesによるファイルパス取得あれこれ。 - Qiita

ディレクトリ配下のファイルの一覧を取得する

読み込ませたいファイルの一覧をリストに格納する

fl <- list.files(path,full.names=T)

※参考にしたサイト

Rでcsvをまとめて読み込む | 分析のおはなし。

読み込み用の関数を作る

例えば、GAのデータであれば、冒頭の6行にコメントが入っているので、それをスキップし、かつCSVファイルの文字コードutf-8なのでそれも指定して読み込む

return (df)を忘れないように

readfiles <- function(x){
  df <- read.csv(x,
                 header=T,
                 fileEncoding = "UTF-8",
                 stringsAsFactors = F,
                 skip=6)
  
  colnames(df) <- c("参照元メディア",
                     "日付",
                    "ページビュー数",
                    "ユーザー")
  df$ユーザー <- as.numeric(sub(",","",df$ユーザー))
  df$ページビュー数 <- as.numeric(sub(",","",df$ページビュー数))
  df$Date <- as.Date(strptime(df$日付,"%Y%m%d"))
  return(df)
}

※その他、Googleアナリティクスのデータはデータが数字にならずにカンマが入っていたり、日付が上手く日付データになっていなかったりのトラップがあるので、読み込み時に対応しておくのが良い。

※参考にしたサイト

qiita.com

関数を使ってファイルをよみこむ

for 分を使って、先程のファイルリストの長さの分だけ繰り返しを適応する。 rbindで繋げていく

df_hoge <- as.data.frame(NULL)
for (i in 1:length(fl)){
  df <- readfiles(fl[i])
  df_hoge <- rbind(df_hoge,df)
}
df <- as.data.frame(NULL)

※参考にしたサイト

blog.okazuki.jp

Rでディレクトリ内のcsvをループで読み込み | バスストップ日誌

Rの時系列分析で、日別の販売データのトレンドを掴む

Excelではどうしても限界が来るからできるだけRを使うべきだというのが持論です。

Rを使えばこんなことも簡易に出来てしまうという例として、周期性のあるデータの可視化があります。

日別販売データは、土日が急に高くなって見辛い

総務省の「家計調査」から

例えば、「家計調査」には日別の集計項目が存在します。

統計局ホームページ/統計FAQ 19A-Q18 1世帯当たり1か月間の日別支出

(総務省統計局トップ > 統計データ > 分野別一覧 > 家計に関する統計 - 家計調査
   家計収支編 > 二人以上の世帯 > 詳細結果表 > 月次 > 表番号6-15 用途分類による日別支出、表番号6-16 品目分類による日別支出)

統計局ホームページ/家計調査(家計収支編) 調査結果

ここから2016年10月-2017年1月の「食料」のデータを引っ張りだして、普通に折れ線グラフにしたらこうなります。

f:id:yyhhyy:20170320193049p:plain

このままでは、いったいいつ頃が増加傾向にあって、いつ頃が減少傾向にあるのか、よくわかりません。

(参考)データを集計しグラフにするまでのコード

dplyr(集計用),reshape2,ggplot2(グラフのため),scales(日付の表示用)などを使う。

youto_10gatsu <- read.csv("youto_10gatsu.csv",header=T)
youto_11gatsu <- read.csv("youto_11gatsu.csv",header=T)
youto_12gatsu <- read.csv("youto_12gatsu.csv",header=T)
youto_1gatsu <- read.csv("youto_1gatsu.csv",header=T)

#データフレームの日付欄をDateデータに変更する関数
datetodate <- function(x){
  x$Date <- as.Date(x$Date)
  return(x)
}

youto_10gatsu <- datetodate(youto_10gatsu)
youto_11gatsu <- datetodate(youto_11gatsu)
youto_12gatsu <- datetodate(youto_12gatsu)
youto_1gatsu <- datetodate(youto_1gatsu)

#データを統合する
youto <- dplyr::full_join(youto_10gatsu,youto_11gatsu)
youto <- dplyr::full_join(youto,youto_12gatsu)
youto <- dplyr::full_join(youto,youto_1gatsu)

#ggplot用にデータを変換
youto_m <- melt(youto,id.vars = c("Date"))


#日付の軸の表示をわかりやすくしグラフにする
g <- ggplot(data=youto_m,aes(x=Date,y=value,color=variable))
g <- g + geom_line(size=1.2)
g <- g + theme_bw()
g <- g + scale_x_date(date_breaks="1 week",
                      date_minor_breaks="1 week",
                      date_labels = "%m/%d")
g <- g + scale_colour_manual(values= c("grey15"))
g <- g + ylab("")
plot(g)

明らかに一週間単位での動きがあるならば、それを除去する作業は簡単

先人の力を借ります。

qiita.com

tjo.hatenablog.com

などを参考にさせて頂きました。

日別データで7日ごとに周期があるのが明白そうであれば、ts関数で季節変動データにする時に「frequency」を「7」にし、stlパッケージでトレンドに分解するとそれらしい分解ができる筈です。

f:id:yyhhyy:20170320194605p:plain

そのコード

#時系列データにする。
youto_t<-ts(as.numeric(youto$食料),frequency=7)
plot(youto_t)

#トレンド等に分解する
youto_t_stl <- stl(youto[,2],s.window="periodic")
plot(youto_t_stl)

せっかくなので、元のデータとトレンドとを一つの線グラフで見てみる

stlパッケージのままでは、x軸が置き換えられてしまって、肝心の何月といった情報が一覧できませんので、改めてトレンド部分を取り出してggplot2などでグラフにするのが良いでしょう。

f:id:yyhhyy:20170320194817p:plain

年末の買い込み以外は平坦ですね。(食料ですから)

今回は地味なオープンデータでしたが、殆どの販売データというものは土日に偏りがあるものですので、日別のデータが手に入る立場にいる人にとっては使う機会がそこそこあるのではないでしょうか?

そのコード

headで適宜データの一部をのぞくと「trend」という項目が、time.seriesの2列めにあることがわかります。それをベクトルデータにして元のデータフレームに列名を付けて追加し、そのデータを元にggplot2でグラフにします。

#トレンドだけ取り出す
head(youto_t_stl)
head(youto_t_stl$time.series)
head(youto_t_stl$time.series[,2])

youto_trend <- youto_t_stl$time.series[,2]

#トレンド項をデータに追加する
youto$trend <- as.vector(youto_trend)

#トレンドをグラフに
youto_t_m <- melt(youto,id.vars =c("Date"))

g <- ggplot(data=youto_t_m,aes(x=Date,y=value,group=variable,color=variable))
g <- g + geom_line(size=1.2)
g <- g + theme_bw()
g <- g + scale_x_date(date_breaks="1 week",
                      date_minor_breaks="1 week",
                      date_labels = "%m/%d")
g <- g + scale_colour_manual(values = c("grey65", "grey15"))
g <- g + ylab("")
plot(g)

何をやっているのか?

7日で周期性のあるとモデル化した部分とそれ以外のトレンド部分とに分けてモデリングしてあてはめているのでしょうか?日本語でググる限り詳解されたサイトはなかなか見つからないので、そう安々と理解できるものではなさそうです。 時系列分析の書籍をあたってアルゴリズムを確認しているうちに日が暮れるので、取り敢えず今は便利なツールとして使っておきましょう。

Rと時系列(3)

uncorrelated.hatenablog.com

Excelにしたければ、

他人と共有する上で、画像のグラフだけだと安心しない人が多いので、csvファイルで吐き出せばExcelで渡すこともできます。

write.csv(youto,"youto.csv",quote=F)

f:id:yyhhyy:20170320195422p:plain

今回のコード全体

library("dplyr", lib.loc="C:/r/R-3.3.2/library")
library("ggplot2", lib.loc="C:/r/R-3.3.2/library")
library("reshape2", lib.loc="C:/r/R-3.3.2/library")
library("scales", lib.loc="C:/r/R-3.3.2/library")
library("RColorBrewer", lib.loc="C:/r/R-3.3.2/library")

youto_10gatsu <- read.csv("youto_10gatsu.csv",header=T)
youto_11gatsu <- read.csv("youto_11gatsu.csv",header=T)
youto_12gatsu <- read.csv("youto_12gatsu.csv",header=T)
youto_1gatsu <- read.csv("youto_1gatsu.csv",header=T)

#データフレームの日付欄をDateデータに変更する
datetodate <- function(x){
  x$Date <- as.Date(x$Date)
  return(x)
}

youto_10gatsu <- datetodate(youto_10gatsu)
youto_11gatsu <- datetodate(youto_11gatsu)
youto_12gatsu <- datetodate(youto_12gatsu)
youto_1gatsu <- datetodate(youto_1gatsu)

#データを統合する
youto <- dplyr::full_join(youto_10gatsu,youto_11gatsu)
youto <- dplyr::full_join(youto,youto_12gatsu)
youto <- dplyr::full_join(youto,youto_1gatsu)

#とりあえずグラフにする
youto_m <- melt(youto,id.vars = c("Date"))

g <- ggplot(data=youto_m,aes(x=Date,y=value,color=variable))
g <- g + geom_line(size=1.2)
g <- g + theme_bw()
g <- g + scale_x_date(date_breaks="1 week",
                      date_minor_breaks="1 week",
                      date_labels = "%m/%d")
g <- g + scale_colour_manual(values= c("grey15"))
g <- g + ylab("")
plot(g)
ggsave(plot=g,file="20170320-1.png",dpi=300,width=4,height=3,scale=2)

#季節データ化
youto_t<-ts(as.numeric(youto$食料),frequency=7)
plot(youto_t)
#分解
youto_t_stl <- stl(youto[,2],s.window="periodic")
plot(youto_t_stl)

#トレンドだけ取り出す
head(youto_t_stl)
head(youto_t_stl$time.series)
head(youto_t_stl$time.series[,2])

youto_trend <- youto_t_stl$time.series[,2]

#トレンド項をデータに追加する
youto$trend <- as.vector(youto_trend)

#トレンドをグラフに
youto_t_m <- melt(youto,id.vars =c("Date"))

g <- ggplot(data=youto_t_m,aes(x=Date,y=value,group=variable,color=variable))
g <- g + geom_line(size=1.2)
g <- g + theme_bw()
g <- g + scale_x_date(date_breaks="1 week",
                      date_minor_breaks="1 week",
                      date_labels = "%m/%d")
g <- g + scale_colour_manual(values = c("grey65", "grey15"))
g <- g + ylab("")
plot(g)
ggsave(plot=g,file="20170320-2.png",dpi=300,width=4,height=3,scale=2)

write.csv(youto,"youto.csv",quote=F)

AI(人工知能)の本領は、「集合知」と「センシング」だと思う

最近はバズワードのように「AI(人工知能)」が流行っている。 とは言え、ニューロ家電やファジー家電と言っていた頃もあるし、一体今、何が重要で、何が重要でないのか?よくわからないというのが正直なところではないかと思う。

人工知能の歴史

人工知能の歴史については、この本が詳しい。

人工知能自体はずっと昔から研究されているもので、これまで何度も流行りがあった。

機械学習、深層学習、人工知能

現在流行っているAIというのは、ドラえもん鉄腕アトムのようなものではなく、「機械学習」「パターン認識」のものである。 感情を理解するお友達ロボットではなくて、与えられた情報を事前に学習した情報と照らし合わせて正しく分類しているだけで、 例えれば、受験勉強が得意な優等生タイプ、ということだ。

機械学習人工知能の言葉の包含関係はこちらの図を見て欲しい。 blogs.nvidia.co.jp

ここ最近発達してきた理由

では、昔からあった筈の人工知能が何故この数年急速に発達してきたのか?というと、その技術的理由は、 - コンピュータの性能が上がって、これまで計算できなかった量の計算ができるようになったこと。 - 乱数によって擬似的に学習用データを自動生成する、という手法を思いついたこと。 である。

人工知能の得意分野を俯瞰で見る

人工知能が人の仕事を奪うのか?誰の仕事を奪うのか?を議論することに意味は薄い。今も昔も技術革新で職を失う人は居たし、これらかも同じだ。 そんな心配をしている暇があったら、人工知能で人間がより豊かになるにはどうしたらいいのか?を考える方が建設的である。

大量のデータが処理できるなら「集合知」が重要である

インターネットで世界が繋がり、Wikipediaが流行りだした頃によく言われていたのが「集合知」という概念である。 Amazonのレコメンド機能が始まったとき、その精度の高さに驚いた頃を思い出して欲しい。過去にその本を買った人が他にどんな本を買っていたのか?という情報を大量に集めれば集める程、「この本を買った人は、この本もチェックしています」の精度が上がる。

となるど、とにかくデータを集めることが重要で、自社だけで頑張って集めるよりも、如何に外部と協力してデータを集めるか?或いは、何かのサービスのついでにデータを収集するか?が人工知能の性能に直結する。

そういう意味で、個人のデータを記録し学習し気の利いたタイミングで何かをサジェストしてくれる人工知能家電、というのは、人工知能のメリットを活かせていないと言えるだろう。

情報が共有されていない分野の情報をつなぐことが重要

今、医療分野に人工知能が入ろうとしている。これまで属人的な個人のスキルとして蓄積されていたものを繋げてみんなで一緒に考えることができるようになる。そういう分野こそ活躍の場があると思われる。

例えば、「小さい子どもの子育て」という分野は、人類が産まれてこの方、うん十年とあるわけだがいつまでたっても正解がなく、毎年のように説は新しくなる。個人にとっては、人生の中で数年間のことで、老人になると必要の薄い知識になる。一方で社会全体では、いつも大量にその情報を欲している人がいる。

そういった問題を解決するために使われるべきだろう。

人間を超えるには「センシング」技術の応用が重要である

このことについて分かり易いのは、落合陽一先生のこの本の最後の方である。

魔法の世紀

魔法の世紀

人間の知覚を超える技術としては、今はおなじみの赤外線が代表的だが、センシング技術が発展することによって、人間の閾値を超えた範囲の空気の振動が測定できるようになった。つまりそれが落合陽一先生の言う「魔法」の世界である。

センシングは既に役立っている

経験と勘に頼っていた農業がセンシングデータによって手軽になってきたことは、周知の事実だと思われる。 もっと身近な例で言えば、デジタルカメラが自動的に人間の顔を認識して自動でピントを併せてくれる、といったものも普通に使っている筈だ。ちょっと前までは人間がマニュアルで焦点距離・光量・シャッタースピード全て調整していた。

人間には感知できないセンシングが重要

相手が怒っているのか? 個人の経験では判断に限界があることがある。

しかし人間よりも細かく声の抑揚を感知できれば、感情も認識できるようになる。

nlab.itmedia.co.jp

「勘のいい人」「気の利く人」すらも人工知能の力でバックアップできる時代が来るかもしれない。

人間がデータとして普段意識していないこともデータベース化する

入館IDの行動ログをとって、誰と誰との会話量が多いか?など組織や個人が普段集計していないデータも収集することによって組織変更に応用できる

或いは、こちらの本にある、バルーンチャレンジの例のように、情報提供のアシストだけをして結果に一見コミットしていないように見えるけれど実際には重要な役割をした人、に注目すると良い場合がある。

ソーシャル物理学:「良いアイデアはいかに広がるか」の新しい科学

ソーシャル物理学:「良いアイデアはいかに広がるか」の新しい科学

こういったことは既に実用化されていて、客先とのメールの回数によって担当営業がそろそろ営業に行くべきクライアントを見つけてアラートを出す、みたいなシステムも存在する。そういった個人では気がつかないデータを活躍することも重要と思われる。

まとめ

  • 色々な人が経験している筈なのに知識が共有されていないなぁという分野
  • 誰かの勘に頼っているなぁという分野
  • データを取ることなんて無理だと思っていた分野

が今後人工知能が伸びる分野だと思われる。

個人の購買履歴を追ってターゲティングされた広告を打つ、或いは、課金ユーザーを維持する、という狭い範囲で統計学が活かされても余り社会を良くするとは思えないし、そういう分野は大きなイノベーションには繋がらない。社会全体の役に立つように統計学を生かして欲しいと思う。