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

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

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

【データ可視化】誤解を招かないように配慮された最新のカラーマップ(カラーパレット)について

先日読んでいた本

科学技術計算のためのPython入門 ――開発基礎、必須ライブラリ、高速化

科学技術計算のためのPython入門 ――開発基礎、必須ライブラリ、高速化

にとても重要なことが記載されていました。

このレインボーカラーマップには問題があるためmatplotlib 2.0からデフォルトのカラーマップがjetからviridisと呼ばれるものに変わります。

「科学技術計算のためのPython」p285

なんと見慣れたあの虹色のカラーマップは今では非推奨なのです。

何が問題なのか

ぶっちゃけわかりにくい。

詳しいことはこちらの記事に説明動画へのリンクがあります。

medvis.org

虹色のjetは確かにカラフルで見栄えこそいいのですが、どの色がどの色より密度が高いことになっているのかないのかピンと来ない、人によってはそれが”ないデータの特徴が存在しているかのように誤解する”ということがあります。

その他にも新しく採用されたviridisは階調の変化が視覚的に分かり易いようにかなり研究されているようです。

色覚異常にも配慮がない

一方で新しくPythonのMatplotlibで採用されたviridisは地味ですが、ちゃんと色の違いがわかることと、所謂、色弱への配慮があります。

色覚異常者がカラーマップを見たときの見え方はこちら。

The viridis color palettes

AndroidもLolipopから色覚異常対応がデフォルトで実装されるようになりました。

androidlover.net

マイノリティへの配慮が出来ないのは、今日日時代に取り残されていると言えるでしょう。

実際にPythonで使ってみる。

色の階調を使って密度を可視化するケースは多々あります。今回は、ヒートマップで使ってみることにしました。

データの取得

総務省の情報通信白書から、時間帯別のメディア行為率を使ってみます。

総務省|平成28年版 情報通信白書|主なメディアの利用時間帯

こちらのcsvファイルを使います。

データの読み込み

Pythonではデータの集計にはpandasが便利です。

引数はこちらのサイトを見て確認しました。

Pandas で CSV ファイルやテキストファイルを読み込む – Python でデータサイエンス

今回のデータは、文字コードシフトJISになっていたので、列名・行名を指定して読み込みます*1

import pandas as pd
heijitsu = pd.read_csv("n5205021_2.csv",encoding="SHIFT-JIS",header=0,index_col=0)

すると、こんな感じのデータフレームです。

テレビ(リアルタイム)視聴 テレビ(録画)視聴 ネット利用 新聞閲読 ラジオ聴取
5時台 6.6 0.4 3.1 2.7 1.0
6時台 22.9 0.4 10.4 7.6 2.0
7時台 31.9 0.5 16.6 8.4 2.4
8時台 20.2 0.5 19.7 5.9 2.9
9時台 8.4 0.9 16.9 2.4 2.2
10時台 6.0 1.2 15.2 1.3 2.2
11時台 6.4 1.2 13.5 0.6 1.9
12時台 14.8 1.1 26.7 1.0 2.0
13時台 9.4 1.9 15.3 0.9 1.7
14時台 6.6 2.2 12.9 0.7 1.7
15時台 6.1 1.5 12.9 0.7 1.2
16時台 7.1 1.4 14.1 0.6 1.3
17時台 12.2 1.9 17.8 2.4 1.7
18時台 19.8 2.0 18.9 1.9 1.6
19時台 34.8 2.3 20.2 2.5 1.1
20時台 41.5 3.0 22.8 3.1 0.9
21時台 42.4 4.8 26.6 2.7 1.0
22時台 35.6 5.0 26.2 1.7 0.6
23時台 18.1 3.6 19.0 0.9 0.5
24時台 6.7 1.8 8.3 0.3 0.4
1時台 2.1 0.7 3.1 0.1 0.4
2時台 0.6 0.2 1.3 0.0 0.2
3時台 0.4 0.1 0.8 0.0 0.2
4時台 0.4 0.1 0.7 0.1 0.3

先ずはデフォルトの色でヒートマップ化

Pythonでグラフを作る時は、seabornを使うと手軽で美しくなります。

qiita.com

今のデータでヒートマップにするには、seabornとmatplotlibからpyplotを呼び出します。また、日本語が豆腐にならないようにフォントを指定します。*2

import seaborn as sns
import matplotlib.pyplot as plt
sns.set(font=['IPAexGothic'])
plt.figure(figsize=(10,10))
sns.heatmap(heijitsu,annot=True)
plt.show()

正直、デフォルトの色も悪くありません。

f:id:yyhhyy:20170504173941p:plain

jetとviridisを比較してみる

試しにjet(raibow)で描いてみます

plt.figure(figsize=(10,10))
sns.heatmap(heijitsu,annot=True,cmap="jet")
plt.show()

赤と黄色とどっちがどれくらい違うのか一発で理解できるでしょうか?

f:id:yyhhyy:20170504174124p:plain

噂の最新のカラーマップviridisを使ってみます。

plt.figure(figsize=(10,10))
sns.heatmap(heijitsu,annot=True,cmap="viridis")
plt.show()

これなら重要なところの方が明るく目立つようになっていると思います。

f:id:yyhhyy:20170504174237p:plain

今回のコード

Anacondaを使っているので「matplotlib inline」を指定しています。

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline

heijitsu = pd.read_csv("n5205021_2.csv",encoding="SHIFT-JIS",header=0,index_col=0)
heijitsu

sns.set(font=['IPAexGothic'])

plt.figure(figsize=(10,10))
sns.heatmap(heijitsu,annot=True)
plt.show()

plt.figure(figsize=(10,10))
sns.heatmap(heijitsu,annot=True,cmap="jet")
plt.show()

plt.figure(figsize=(10,10))
sns.heatmap(heijitsu,annot=True,cmap="viridis")
plt.show()

yyhhyy.hatenablog.com

yyhhyy.hatenablog.com

Rにもあるらしい

ほんとこのサイトは情報が早いですね。

www.karada-good.net

*1:元のCSVから事前に何行か削除しています。

*2:Macの場合はOsakaなどが無難でしょう

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

広告業界ではある媒体に出稿した時の接触回数の分布(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)