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

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

アンケートや視聴率の誤差を推測するプログラム(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の行動ログをとって、誰と誰との会話量が多いか?など組織や個人が普段集計していないデータも収集することによって組織変更に応用できる

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

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

まとめ

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

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

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

当選確率nのくじをn回ひいた時に少なくとも1回は当る確率の計算について

先日このツイートを見て、念のためそれを確かめてみた。

n回中、少なくとも1回は当る確率

「n回の試行中、少なくとも1回は当たりが出る確率」は、「n回の試行中、1回も当たらない確率」を「1(=100%)」から引けばいい。 *1

当る確率が「p」だとして、1回の試行で当たらない確率は「1- 1/p」。 n回連続で当たらない確率は、「(1-1/p)n」(累乗計算)となります。

試しにPythonで数式計算をしてみる

Pythonで数式計算をするライブラリは「sympy」

from sympy import *

sympyでは先ず使用する変数を定義します。 ここでは1回の当選確率を「p」。求めたい値であるn回試行中1回は当る確率を「p1」。試行回数を「n」としました。

n,p,p1=symbols("n p p1")

先ほどの計算式を記述

p1 = 1 - (1-1/p)**n

また、今の設定では、当選確率の逆数の回数分だけくじを引くことになっているので「p == n」を代入し、それを新たに「p1_2」と定義します。

p1_2 = p1.subs(p,n)
print(p1_2)

すると、「1/n」の確率で当るくじを「n」回ひいた時に、少なくとも1回は当るが出る確率は

「-(1 - 1/n)**n + 1」

とわかりました。 先ほどの手計算の数式と同じです。

このp1_2に具体的に「n==2」なり「n==100」なりを代入した結果を出力すれば、求めたい数値が出てきます。 「evalf()」で分数が小数点表示になります

print(p1_2.subs(n,2).evalf())
print(p1_2.subs(n,6).evalf())
print(p1_2.subs(n,100).evalf())

当選確率1/2で2回だと「0.750000000000000」 1/6の確率で6回だと「0.665102023319616」 1/100の確率で100回だと「0.633967658726771」

値を求めるだけならRでできる

数式を求めることに意義が薄いと思えば最初からRを使えば手っ取り早いです。

当たり/外れ の2種類しかない試行は二項分布になります。 確率を求めるのは、「d~」で、二項分布は「binom」。 引数は、「当たりの確率」、「試行回数」、「当る確率」の3つで、先ほど同様、1回も当たらない(==0)確率を1から引けば良いです。

1-dbinom(0,2,1/2)
1-dbinom(0,4,1/4)
1-dbinom(0,6,1/6)
1-dbinom(0,100,1/100)
> 1-dbinom(0,2,1/2)
[1] 0.75
> 1-dbinom(0,4,1/4)
[1] 0.6835938
> 1-dbinom(0,6,1/6)
[1] 0.665102
> 1-dbinom(0,100,1/100)
[1] 0.6339677

それはそうと、さっきと数字が違うじゃないですか?

さっきのツイートの数字は、63%と断言されていましたが、実際には違う、、、でも1/100で100回の時はとても近い

Rで繰り返し処理

というわけで、n==2 から n==100 までの夫々の確率をデータフレームに格納して、プロットしてみます。

繰り返し処理でデータフレームに追記していく場合、rbindを使います。

df <- data.frame()
for (n in 2:100){
  dat <- data.frame(1- dbinom(0,n,1/n))
  df <- rbind(df,dat)
}
colnames(df) <- c("prob")
df$x <- rep(2:100)

結果を出力すると、試行回数が上がるに連れて一定の数字に修練しているようです。

> head(df)
       prob x
1 0.7500000 2
2 0.7037037 3
3 0.6835938 4
4 0.6723200 5
5 0.6651020 6
6 0.6600833 7
> tail(df)
        prob   x
94 0.6340653  95
95 0.6340450  96
96 0.6340250  97
97 0.6340055  98
98 0.6339864  99
99 0.6339677 100

グラフにしてみる

ggpotでグラフ化

p <- ggplot(data=df,aes(x=x,y=prob))
p <- p + geom_bar(stat="identity")
p <- p + ggtitle("当たりの確率の逆数と同じ回数だけひいたら1回は辺りが出る確率")
p <- p + geom_text(aes(x=100,
                       y=df[df$x==100,]$prob,
                       label=round(df[df$x==100,]$prob*100,2)),
                   vjust=-1)
plot(p)
ggsave(plot=p,file="20170115.png",h=4,w=6,scale=2)

ほぼほぼ63.4%に収斂しています。

f:id:yyhhyy:20170115172459p:plain

それでは、「3n回」や「5n回」は?

繰り返し計算中の確率計算部分で試行回数を「n3」や「n5」にするだけです

dat <- data.frame(1- dbinom(0,n*3,1/n))

3倍の回数なら

f:id:yyhhyy:20170115172557p:plain

5倍の回数なら

f:id:yyhhyy:20170115172609p:plain

実際のところ、確率1%であっても100回も引く人はいないと思いますが、、、

これだけ回数を増やせばポアソン分布である

ガチャの確率が一般的にどれくらいなのか知りませんが、 二項分布は確率が低く、試行回数が増えればポアソン分布で近似できる、ということがわかっています。

二項分布 - Wikipedia

dskjal.com

qiita.com

例えば、当選確率pが「1/100」で試行回数nが「100回」なのであれば、

lambda = n * p は、「100*(1/100)」

同様に「300回」なら「100(1/300)」、「500回」なら「500(1/100)」を引数として

1 - dpois(0,(100*(1/100)))
1 - dpois(0,(300*(1/100)))
1 - dpois(0,(500*(1/100)))

とすれば、先ほど収束していった値が出てきます。

> 1 - dpois(0,(100*(1/100)))
[1] 0.6321206
> 1 - dpois(0,(300*(1/100)))
[1] 0.9502129
> 1 - dpois(0,(500*(1/100)))
[1] 0.9932621

コード全体

from sympy import *
n,p,p1=symbols("n p p1")
p1 = 1 - (1-1/p)**n
p1_2 = p1.subs(p,n)
print(p1_2)
print(p1_ans.evalf())
print(p1_2.subs(n,6).evalf())
print(p1_2.subs(n,100).evalf())
library("dplyr", lib.loc="C:/hogehoge/r/R-3.3.2/library")
library("reshape2", lib.loc="C:/hogehoge/r/R-3.3.2/library")
library("ggplot2", lib.loc="C:/hogehoge/r/R-3.3.2/library")

1-dbinom(0,2,1/2)
1-dbinom(0,4,1/4)
1-dbinom(0,6,1/6)
1-dbinom(0,100,1/100)

df <- data.frame()
for (n in 2:100){
  dat <- data.frame(1- dbinom(0,n,1/n))
  df <- rbind(df,dat)
}
colnames(df) <- c("prob")
df$x <- rep(2:100)
head(df)
tail(df)

p <- ggplot(data=df,aes(x=x,y=prob))
p <- p + geom_bar(stat="identity")
p <- p + ggtitle("当たりの確率の逆数と同じ回数だけひいたら1回は辺りが出る確率")
p <- p + geom_text(aes(x=100,
                       y=df[df$x==100,]$prob,
                       label=round(df[df$x==100,]$prob*100,2)),
                   vjust=-1)
plot(p)
ggsave(plot=p,file="20170115.png",h=4,w=6,scale=1)

df_3 <- data.frame()
for (n in 2:100){
  dat <- data.frame(1- dbinom(0,n*3,1/n))
  df_3 <- rbind(df_3,dat)
}
colnames(df_3) <- c("prob")
df_3$x <- rep(2:100)
print(df_3)

p3 <- ggplot(data=df_3,aes(x=x,y=prob))
p3 <- p3 + geom_bar(stat="identity")
p3 <- p3 + ggtitle("当たりの確率の逆数との3倍の回数をひいたら1回は辺りが出る確率")
p3 <- p3 + geom_text(aes(x=100,
                       y=df_3[df_3$x==100,]$prob,
                       label=round(df_3[df_3$x==100,]$prob*100,2)),
                   vjust=-1)
plot(p3)
ggsave(plot=p3,file="20170115_3.png",h=4,w=6,scale=1)


df_5 <- data.frame()
for (n in 2:100){
  dat <- data.frame(1- dbinom(0,n*5,1/n))
  df_5 <- rbind(df_5,dat)
}
colnames(df_5) <- c("prob")
df_5$x <- rep(2:100)
print(df_5)

p5 <- ggplot(data=df_5,aes(x=x,y=prob))
p5 <- p5 + geom_bar(stat="identity")
p5 <- p5 + ggtitle("当たりの確率の逆数との3倍の回数をひいたら1回は辺りが出る確率")
p5 <- p5 + geom_text(aes(x=100,
                         y=df_5[df_5$x==100,]$prob,
                         label=round(df_5[df_5$x==100,]$prob*100,2)),
                     vjust=-0.5)
plot(p5)
ggsave(plot=p5,file="20170115_5.png",h=4,w=6,scale=1)

1 - dpois(0,(100*(1/100)))
1 - dpois(0,(300*(1/100)))
1 - dpois(0,(500*(1/100)))

*1:数学が面白いのはこの手の逆転の発想が普通に出てくることですよね

【夢を見た】VMWareとEl Capitanのメモ

夢を見た話。

以前に、

yyhhyy.hatenablog.com

という記事を書きましたが、今はVMwarePlayerを使用しています。

PCを新しくしたのでこれを機会にMavericksからEl Capitanに変更したので、その時に参考にしたリンク集。これらに従えば簡単に使用できるようになります。

参考サイト一覧

El CapitanのISOイメージファイルを手に入れるまで

MacOSはISOイメージの形で頒布されていないので以下の手法でISOイメージファイルを作る

Windows で VMWare に Mac OSX El Capitan を入れるには - Enjoi Blog

※当然ですがMacOSを持っていないとMaxOSはダウンロードできません。

その際の画面等

基本はソフトウェアの指示の通り。

seedofhack.wp.xdomain.jp

MacOSインストール後

VMWare toolsのインストール方法

英語ですがここがわかりやすいです。

www.wikigain.com

キーボードの配列の有効化

(*何故か私は上手く行きませんでした。)

ichitaso.com

LinuxVMware

インストールまで

そもそもLinuxVMwareを入れる方法がわかりにくい

itlx.ldblog.jp

itlx.ldblog.jp

インストール後

バーチャルウィンドウは、Ctrl + Alt + Enter で全画面表示に切り替わりますが、 微妙に画面の端っこをクリックするとUbuntuのLauncherの方が反応してしまいます。

UbuntuのLauncherは、このように下側に配置し、

Ubuntu 16.04 LTSのLauncherを画面下に配置する

MacのDockの方は、Dockの環境設定で右か左にしておけば、操作中に誤ってホストOS側に移ってしまう、ということもないでしょう。