【備忘録】_PythonでGoogleトレンドの季節変動を取り除いてみる
自分の備忘録目的ので、そりゃそうだよね？ということがただ確認できる、というもので、特に発見を目的としたエントリーではありません。
Googleトレンドから「ドラえもん」のデータを取得
Googleトレンドでエリアを日本に絞り、過去5年の「ドラえもん」の検索数を取得します
Jupyterで季節変動をとる
データの読み込みからプロットまで
こちらを参考に致しました。 data.gunosy.io
必要なパッケージの読み込み。
import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns %matplotlib inline
Googleトレンドでは3行目に列名がありますので、頭2行はスキップして読み込みます
df = pd.read_csv("doraemon.csv",skiprows=2)
renameで列名を英語にします。*1
df = df.rename(columns = {"週":"week","ドラえもん: (日本)":"doraemon"})
脚書は日付は文字列として認識されているので日付データ（dattime）に変換します
df["week"] = pd.to_datetime(df["week"])
| week | doraemon | |
|---|---|---|
| 0 | 2012-06-17 | 21 |
| 1 | 2012-06-24 | 20 |
| 2 | 2012-07-01 | 21 |
| 3 | 2012-07-08 | 22 |
| 4 | 2012-07-15 | 23 |
グラフを確認すると緩やかな山型のようです。
plt = df.plot()
月別に集計し直す
Googleトレンドは週単位ですので、これを月あたりの合計値に再集計します。*2
groupby でできる、resample でできるなど、色々と情報はありましたが、上手くいかず、このサイトに従いました。
df_m = pd.pivot_table(df,index=pd.Grouper(key="week",freq="M"),values="doraemon",aggfunc=sum)
week
2012-06-30 41
2012-07-31 111
2012-08-31 103
2012-09-30 188
2012-10-31 86
Name: doraemon, dtype: int64
pandasのデータフレームに戻します。
df_m = pd.DataFrame(df_m)
| doraemon | |
|---|---|
| week | |
| 2012-06-30 | 41 |
| 2012-07-31 | 111 |
| 2012-08-31 | 103 |
| 2012-09-30 | 188 |
| 2012-10-31 | 86 |
図示します
df_m.plot()
先程と余り印象は変りません。
月単位のトレンドを取り除きます。
1つ目のブログに従います。satsmodelsというライブラリを使います。
import statsmodels.api as sm
12ヶ月周期があると仮定することで月の変動を取り除けます。
df_ts = sm.tsa.seasonal_decompose(df_m.values,freq=12)
df_ts.plot()
上が元々のデータ。残差が’一番下ですが、残差もそこそこあるようです。。。
季節変動を除いたトレンドを確認
トレンドは「trend」に入っているので、先程のデータの日付のインデックスと併せてデータフレームにします。
df_trend = pd.DataFrame({"trend":df_ts.trend,"date":df_m.index})
df_trend = df_trend.set_index(["date"])
| trend | |
|---|---|
| date | |
| 2012-06-30 | NaN |
| 2012-07-31 | NaN |
| 2012-08-31 | NaN |
| 2012-09-30 | NaN |
| 2012-10-31 | NaN |
個別に図示するとこうなります。
plt= df_trend.plot()
2014年～2015年の盛り上がりはなんだろう？と思い返すと、あの時、「STAND BY ME ドラえもん」という映画がヒットしましたね。あと、この頃からドラえもんの人形がプロモーションで東京タワーやらなんやらに出没するようになった記憶があります。
月別の変動を確認
同様に今度はトレンドではなく季節変動の方を取り出します
df_seasonal = pd.DataFrame({"seasonal":df_ts.seasonal,"date":df_m.index})
df_seasonal = df_seasonal.set_index(["date"])
| seasonal | |
|---|---|
| date | |
| 2012-06-30 | -14.965972 |
| 2012-07-31 | -10.278472 |
| 2012-08-31 | 57.315278 |
| 2012-09-30 | -5.538889 |
| 2012-10-31 | -29.122222 |
グラフより表の方がみやすいかもしれません。ドラえもんお盛り上がりは、8月と3月のようです。3月は毎年新作映画が公開されますね。
夏って、何でしたっけ？
というわけで落ちもなにもありませんが、以前、RでやったことをPythonでもやってみました。
csvファイルで結果を保存すれば、これで長期的なGoogleトレンドによるトレンドを抜き出しました。
df_trend.to_csv("doraemon_trend.csv")
時系列データはマーケティングでは頻繁に出てくると思うのですが、余り聞かないのは、実際のところ時系列データをそのままモデリングするのは結構たいへんだからだと思います。
とはいえ、ざっくりとしたトレンドと季節変動は分解して図示した方が、頭の整理はされそうな気がします。幸い作業自体は楽なので。
最近AKBのミュージックビデオの出来がとてもいいと思った
たまたまYouTubeで他の曲を聞いていたらAKBの「ハイテンション」のPVが”次の動画”に出てきて、聞いてみたら（というか見てみたら）何か凄くよかった。
何がそんなによかったのか？整理できていなけれど、こういうファーストインプレッションは、後で当たり前になってしまい忘れてしまいがちなので、取り敢えず書きなぐっておきたい。
「ハイテンション」について
- 出だしのパッチンでライトがつくのに気がつく、という入り方がファンタジックでいい。（オチで現実に戻っているも洒落ている）
- パッチンをキッカケに画面が切り替わってモブが踊り出す緩急がいい
- 服のカラーコーディネートがいい。ちょっと昔のモダンガールみたいな色とか、目立つ色がいい。
- ちょっと俯瞰からのカメラでダンスを撮ってから個人のアップを抜くのがいい
- 人差し指を振る謎の振り付けがジワジワカワイイ
- 古典的だがエンドロールで楽しそうな撮影風景が入るのもいい
その他のMVについて
歌っている人が笑顔のPVは見ていて気持ちがいい
ハイテンションもそうだけど、恋するフォーチュンクッキーのときもそうだった。
【MV full】 恋するフォーチュンクッキー / AKB48[公式]
真面目な顔で歌っているMVもいいけど、楽しそうに歌っている姿というのは、それだけでポイントがあがるものだと思う。
この曲は、上2曲ほど分かりやすい笑顔ではないけれど、凄く楽しそうにしている笑顔がちょいちょい挟まって、全体的に伸び伸びと歌っている印象が伝わる。
【MV】Set me free Short ver.[TeamA] / AKB48[公式]
ノスタルジックな演出がいい
メロディはどっかで聴いた感じもあるけど、この女学校やスケバン時代のノスタルジックなモチーフがいい。
【MV】アクシデント中 Short ver.〈AKB48 U-19選抜〉/ AKB48[公式]
「逆さ坂」の歌謡曲からスタートするのも味がある。
【MV】逆さ坂 Short ver.〈じゃんけん民〉/ AKB48[公式]
哀愁のトランペッターは歌詞がノスタルジックなのだが、歌われる情景が映像に思い浮かぶ。（キャバレーに居る感じしません？）
【MV】哀愁のトランペッター Short ver.[TeamK] / AKB48[公式]
余談だが、花の垣根から顔だけ出しているのも面白い。
アイドルというより大人の色気が映像で映えるようになってきた気がする
最初の頃はよく水着になったり、そんな映像がYouTubeで再生回数が伸びるのなんか狡いだろー！と思っていたけど、最近はそうでもない渋いMVが増えてきた。
特に、山本彩と、この数年でちょっと痩せた渡辺麻友は、映像映えがする。
オジサンキラー感がある
学校を卒業して戻ってくる演出とか、ノスタルジックとは別の次元でオジサンキラーなテーマの選び方が多い。
【MV】光と影の日々 Short ver. / AKB48[公式]
おばあちゃんから貰った500円玉を使わないで持っていた、とか、そういうテーマも凄くオジサンキラーである。そう言えば、そんな感覚、あの頃はあったよねと。秋元康はこういうところが天才的だと思う。
【MV】あの頃の五百円玉 Short ver. / AKB48[公式]
ドキュメンタリータッチの演出がベストマッチするようになってきた
これはアイドルグループのMVだろうか？？当然、取材先の同意があるかるからできるのだと思うので、いつの間にか地方の学校の人もOKを出すような非アンダーグラウンド感がAKB48にはあるということだ。
【MV full】願いごとの持ち腐れ / AKB48[公式]
なかなかまとまらないけれど、他の数あるアイドルグループとは少し違う感じのミュージックビデオなのだ。
【データ可視化】誤解を招かないように配慮された最新のカラーマップ（カラーパレット）について
先日読んでいた本
![科学技術計算のためのPython入門 開発基礎、必須ライブラリ、高速化 [ 中久喜健司 ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/3886/9784774183886.jpg?_ex=128x128 "科学技術計算のためのPython入門 開発基礎、必須ライブラリ、高速化 [ 中久喜健司 ]")
にとても重要なことが記載されていました。
このレインボーカラーマップには問題があるためmatplotlib 2.0からデフォルトのカラーマップがjetからviridisと呼ばれるものに変わります。
「科学技術計算のためのPython」p285
なんと見慣れたあの虹色のカラーマップは今では非推奨なのです。
何が問題なのか
ぶっちゃけわかりにくい。
詳しいことはこちらの記事に説明動画へのリンクがあります。
虹色のjetは確かにカラフルで見栄えこそいいのですが、どの色がどの色より密度が高いことになっているのかないのかピンと来ない、人によってはそれが”ないデータの特徴が存在しているかのように誤解する”ということがあります。
その他にも新しく採用されたviridisは階調の変化が視覚的に分かり易いようにかなり研究されているようです。
色覚異常にも配慮がない
一方で新しくPythonのMatplotlibで採用されたviridisは地味ですが、ちゃんと色の違いがわかることと、所謂、色弱への配慮があります。
色覚異常者がカラーマップを見たときの見え方はこちら。
AndroidもLolipopから色覚異常対応がデフォルトで実装されるようになりました。
マイノリティへの配慮が出来ないのは、今日日時代に取り残されていると言えるでしょう。
実際に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を使うと手軽で美しくなります。
今のデータでヒートマップにするには、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()
正直、デフォルトの色も悪くありません。
jetとviridisを比較してみる
試しにjet（raibow）で描いてみます
plt.figure(figsize=(10,10)) sns.heatmap(heijitsu,annot=True,cmap="jet") plt.show()
赤と黄色とどっちがどれくらい違うのか一発で理解できるでしょうか？
噂の最新のカラーマップviridisを使ってみます。
plt.figure(figsize=(10,10)) sns.heatmap(heijitsu,annot=True,cmap="viridis") plt.show()
これなら重要なところの方が明るく目立つようになっていると思います。
今回のコード
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()
Rにもあるらしい
ほんとこのサイトは情報が早いですね。
バナー広告の接触回数分布を関数で回帰してみる
広告業界ではある媒体に出稿した時の接触回数の分布（1回当たった人がn人、2回当たった人がn人）をベータ二項分布などで近似しシミュレーションするのが一般的と言われている
なおこのサイトによるとWeb広告の場合は、負の二項分布がグローバルで標準らしい。*1
ということで、バナー広告の接触回数分布がどんな関数になるのか、入手可能なデータで確かめてみた。
ビデオリサーチ社が2年前に調査したデータを使う
データソースについて
説明はここを参照して欲しい。
リリースの最後の方に、接触回数の分布がある。
オンライン広告の認知効果の基準値整備に関するお知らせ | ニュース || コーポレートサイト
このサンプル数がわからなかったが、計算の簡単のためこのパーセンテージの小数点が消えるように桁を変えて,仮にこの値を回帰できる関数を探してみる。
x y 1 271 2 179 3 129 4 89 5 69 6 53 7 40 8 32 9 25 10 20
プロットしてみる
プロットした時の関数は省略する。
一般的に言われている負の二項分布っぽい感じはする。
負の二項分布について
負の二項分布の理解
負の二項分布は、確率pの試行をN回繰り返す時に、r回成功するまでに何回試行を繰り返すか？の分布を知る為に使われるのだが、今回の場合はその理解というよりも、こちらのブログの通り、
http://sugisugirrr.hatenablog.com/entry/2016/01/31/%E7%AC%AC%EF%BC%97%EF%BC%94%E5%9B%9E_%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E5%88%86%E5%B8%83%E3%81%A8%E3%82%B5%E3%83%83%E3%82%AB%E3%83%BC%E3%81%AE%E3%82%B4%E3%83%BC%E3%83%AB%E6%95%B0sugisugirrr.hatenablog.com
カウントデータの回帰にポワソン分布を使おうとしたら分散がもっと大きい分散なので、ガンマ分布とポアソン分布の階層モデルを作るとき、それが負の二項分布になる、というポアソン分布の過分散時の推定、と言われている方の用途で使う。
このサイトを見ると確かにガンマ分布で何パターンか作ったポアソン分布を合成していくと負の二項分布になっている。。。
http://stat.biopapyrus.net/probability/negative-binomial.html
RStanに負の二項分布のあてはめ
Rのパッケージ、nbGLMでやってもほぼおなじ結果なのだがせっかくなのでRStanでベイズモデリングしてみる。
知識不足のため、収束させる為の条件はわからなかったのでコードはほぼほぼこれらのサイトを見よう見まねで使った。
Simulate and fit negative binomial GLMs in 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を使って値を代入すればよい。
![実践ベイズモデリング 解析技法と認知モデル [ 豊田秀樹 ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/2206/9784254122206.jpg?_ex=128x128 "実践ベイズモデリング 解析技法と認知モデル [ 豊田秀樹 ]")
のコードを見て、”なるほど”と思った。流石です。
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を使う人ととそうでない人とバラつきがあるということだろう。
媒体数が限られているマス広告と比べて、果してバナー広告の接触回数分布が簡単に回帰できる類いのものなのか？という疑問はあるし、そもそもターゲティングが違えばまた違うだろうし、そもそもリマーケティングを入れたりフリークエンシーキャップをかけたらもっと複雑になるので予測すること自体が難しいのではないか？とは思う。
ただ、「負の二項分布で近似できる」という話については、嘘ではなさそうだ。
アンケートや視聴率の誤差を推測するプログラム（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
参考にしたサイト
自分のPCはプロキシの問題があってエラーが出てしまったのでAnacondaをPython2.7のAnaconda2にまで下げて以下を実行しています。
pyinstallerを手に入れる
スタンドアローン化するパッケージは幾つかあるようだがpyinstallerが一番手軽。
以下のサイトの夫々のページのアドバイスに沿ってpyinsatallerを導入する
参考にしたサイト
PythonでWindowsアプリケーションを作ってみよう – Nobwak's Lair
Windowsにpyinstallerをインストール – Nobwak's Lair
pywin32が必要
先程のリンクの２つめ。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を使う。
参考にしたサイト
コマンド自体はとても簡単。
pyinstaller hogehoge.py --onefile --clean
- 単独のファイルにするには「–onefile」のオプションを。
- やり直しする場合は、「–clean」で前のファイルを必ず消す。
- コマンドプロンプトを立ち上げないようにするには、「–noconsole」をつける
なお、実行には10分位は待つ。
そして実行ファイルをダブルクリックしてからかなり待つ
とても便利だけど普通の人にはちょっと待ちすぎ？という位長い時間起動にかかるのがたまにきず。ファイルも重い。
参考にしたサイト
www.slideshare.net
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)
とすればよい
参考にしたサイト
アンケートや視聴率の誤差を推測するプログラム（2）-プログラムのGUI化-
このエントリーは3回中の2回目です
Pythonが入っているPCであれば先般の記事のプログラムで充分なのですが、コマンドプロンプトの画面に数字を打ち込むのは一般の人には抵抗があると思われます。
余談ですがコマンドプロンプトの背景は白の方が心が落ち着くのでオススメです！
というわけで、GUIベースのプログラムに書きかえてみます。
Tkinterを使う
基本的な使い方はこちらの通り
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」が必要です
参考にしたサイト
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)
参考にしたサイト
関数を適用させする
クリックで実行するように、ボタンに対して「command」で先程の関数を適用させます。
main_Button = Tk.Button(text = u"入力したらこのボタンをクリックして下さい",command = kakunin)
main_Button.grid()
参考にしたサイト
コード全体
順番にウィジェットを並べて配置していき、最後に表示させます。
# 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()
数値を入力してクリックするとこのような画面で表示されます
アンケートや視聴率の誤差を推測するプログラム（1）-標本誤差の計算まで-
アンケート調査や視聴率の数字を見せた時に、 「これはn数が少ないから信憑性が薄いねぇ」 と言われることがあると思います。しかし実際のところ、 ”では、今のn数だとどれくらいの幅で理解すべきなのですか？” という質問に応えてくれる人は少ないでしょう。というのも意外とその場で計算するのは容易ではないからです。
今回は、サンプルサイズから標本誤差を推定するスタンドアローンな実行ファイルを作るまで、3回に分けてブログにします。
標本誤差について
どんな統計の本にもあるような基礎的な話ですが、国勢調査などを除いて、殆どの調査は全数調査はしません。
サンプルを集めて、該当質問の平均値を取り、それを結果の値とします。
更に、以下にランダムにサンプルが集められていたとしても偏りが発生していないという絶対的な保証はありません。そのような時、 少なくとも平均値については、サンプル数が多ければ大数の法則がなりたつ ということを使って、 同じような調査を何度も何度もやったらそのうち95％は平均的に○○～○○の値に収まるであろう という何とも消極的なことしか分からないのです。
標本誤差とは？
どういう計算をするのか？はこのサイトが分かりやすいでしょう。
ビデオリサーチ社の視聴率の場合
サンプルサイズから結果の値の標準誤差を知る為のプログラムを作る
- 普段統計の計算はRを使っていますが、今回は最終的に実行ファイルにしたいのでPythonにします。
- 現時点では実行ファイルにするにはPython2の方が良いのでPython2を使います
Pythonを使った標本誤差の計算
正しい手法はこのサイトの通り
「Pythonで正規分布の平均値の信頼区間を計算する方法 」
しかし先程の参照サイトの通り、母数に対してサンプルサイズがすくない場合はt分布ではなく正規分布でも良い、とありましたので、簡略化することにします。
Pythonのscipy statsを使う
使うモジュールは、scipyの中のstatsから、norm.intervalという関数を使って、平均値0、分散1の正規分布の時の95％の有意点の上下を求めます
参考にしたサイト
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)
サンプルサイズやアンケート結果の値を画面から入力する
サンプルサイズをコード本文にベタ打ちすると、毎度毎度コードの修正が必要になってしまうので、コンソールから都度都度入力するようにしましょう
参考にしたサイト
print(u"サンプル数を入力して下さい") n_samples = float(input('>>> ')) print (u"パーセンテージを入力して下さい。（5.5パーセントなら、5.5と入力）") p_ori = float(input('>>> '))
このような感じで自分で数字を打込み、それが引数となるようになります
誤差の値を計算する
一番最初の参照サイトの簡略化された数式を実行します。配列のためだけに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"単位パーセント")
