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

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R: シニングによる非定常ポアソン過程のシミュレーション

ポアソン過程は適当に点を間引いてやっても(非定常)ポアソン過程になる。この性質を使って非定常ポアソン過程をシミュレートすることができて、それをシニング(thinning)と呼ぶ。詳しくは An Introduction to Statistical Computing: A Simulation-based…

バスモデル:エクセルによる需要予測

バスモデルを使って需要予測ができるエクセルシートを以下に公開しました。http://zisatsu.web.fc2.com/monooki/excel/Bass_model.xlsx 使い方 データの列に予測したい対象のこれまでの販売数を入れます。今回は例として トヨタ自動車、ハイブリッド車のグロ…

{googleAnalyticsR} {ggplot2} 何秒でユーザが離脱したか、残存率の図示

直帰率や平均滞在時間よりも詳しくエンゲージメントを把握したい場合、次の図がおすすめです。何秒目に何パーセントのセッションが継続しているかを知ることができます。ここでは新規訪問とリピーター訪問で層別していますが、ランディングページや流入経路…

{googleAnalyticsR} {ggplot2} Crossbar plot で滞在時間の分布を比べる

Crossbar plot は次のような図です。これがなにを表すかは場合によります。今回は箱の上の辺にデータの75%点、真ん中の線に中央値、箱の下の辺に25%点を割り当てました。いわば退化した箱ひげ図です。次の図はあるランディングページを経由した訪問の、月ご…

{ggplot2} 集計済みデータのヴァイオリンプロット

データがこのようにすでに集計された状態で与えられたときに、複数カテゴリ間で分布を比較したい。 value A B 1 1 1 3 2 2 2 2 3 3 3 1 ヴァイオリンプロットが使えることがわかった。 library(cowplot) library(tidyr) data <- data.frame(value = 1:3, A =…

{googleAnalyticsR} {ggplot2} 内訳をドリルダウンしていく棒グラフ

ggplot2 で左から右に内訳をドリルダウンしていく棒グラフ(ツリーマップ?)を書きました。新規訪問で自然検索経由のデスクトップを利用している男性の25歳〜34歳のCVが多いことがわかります。ついで再訪問でソーシャル経由訪問のモバイルを利用している男…

R の treemapify パッケージで Google アナリティクスのデータをツリーマップに

ただパッケージ動かしてみたっていうだけなんですけど。インストール方法は GitHub - wilkox/treemapify: ggplot2 geoms for drawing treemaps に従ってください。ggplot の記法が使えて、facet とかできる点がいいと思いました。 #library(devtools) #inst…

Excel でワイブルプロット

ワイブルプロットを書くためのエクセルファイルを以下に公開しました。http://zisatsu.web.fc2.com/monooki/excel/WeibullPlot.xlsxいまどきだれが使うんだ。ワイブルプロットのなんたるかについてはウィキペディアなどを参照してください。ワイブル分布 - W…

検定いらずの AB テスト:φ 係数を用いたサンプルサイズ設計

AB テストと呼ばれるものや、AB テストをやる目的は様々だろう。ここではウェブ系施策において、「A と B を比べてみて、どっちか良さそうなほうを使いたい」という目的で行われる AB テストで、かつ CTR(クリック率)とか CVR(コンバージョンレート)とか…

[RStan]差分方程式で呂布カルマのフリースタイルダンジョン出場を振り返る

モチベーション [dlm]状態空間モデルでトレンドと広告の効果を分離して推定する - 廿TT で、「広告の効果測定において、残存効果、タイムラグをモデルに組み込みたい」というコメントを頂戴したので、それっぽいモデルを考えてみた。提案モデルはちょっと…

クラッシュロワイヤル進捗

はじめてジャングルアリーナまでいきました。 主に使っているデッキはこれです。

[RStan]同時確率に基づく検索キーワードのクラスタリング(改訂版)

[RStan]同時確率に基づく検索キーワードのクラスタリング - 廿TT の改訂版です。上エントリはそもそもやろうとしてることに無理があったと思い、ある程度クラスタを人力で与えることにしました。 モデル SEO 的には本当は語順にも意味があるんだろうけれど…

[RStan]同時確率に基づく検索キーワードのクラスタリング(失敗版)

このモデルは推定するたびに結果が変わることが判明しました。申しわけありません。以下の記述はなかったことにしてください。改訂版を書きました → [RStan]同時確率に基づく検索キーワードのクラスタリング 2 - 廿TT はじめに 検索キーワードのグルーピン…

R でシミュレーションする格子の上の SIR モデル

SIR モデルは感染症の流行過程のモデルで、詳しくは 感染症のモデル(SIRモデル)に入門した - 廿TT に書いた。あたらしい人工知能の教科書 プロダクト/サービス開発に必要な基礎知識作者: 多田智史,石井一夫出版社/メーカー: 翔泳社発売日: 2016/12/17メデ…

[RStan]項目反応理論の応用でフリースタイルダンジョン登場ラッパーの強さをランキングしてみた

分析対象 フリースタイルダンジョンはフリースタイル(即興)のラップバトルで、チャレンジャーがモンスターと呼ばれる強豪ラッパーを勝ち抜き、賞金獲得することを目指すテレビ番組です。データは、 フリースタイルダンジョン 結果 勝敗 全試合 - 戯言 から…

Google アナリティクスとサーチコンソールのデータをマージして使う

統計ソフト R から Google アナリティクスのデータを引っぱってこれるパッケージ googleAnalyticsR の使い方は、 googleAnalyticsR の使い方(Version:0.1.0) - 廿TT にざっと書いた。サーチコンソールのデータを引っぱってこれるパッケージ searchConsoleR…

Albert (2008) 打者の調子の波のモデル化(後編)

Albert (2008) 打者の調子の波のモデル化(前編) - 廿TT の続きです。以降の分析の目的は、よく言われる野球選手の「調子の波」を選手間で比較可能な指標にするにはどうするか、ということです。分析対象のデータはカルロス・ギーエンという選手の2005年の…

SIR モデルと非定常ポアソン過程

疫学の SIR 動態モデルへのいろいろなあてはめ…MCMCとか粒子フィルターとか Approximate Bayesian Computation (ABC) とか教えてくれるコース https://t.co/vD5Fk9goGQ があるそうで…昔は最小二乗法でも使ってたのかな? https://t.co/tglcj91WtJ— 久保拓弥 …

紀元一世紀のできごとを覚えるために語呂合わせを考えてみました(1〜30)

西暦1年 マヤで供物(人身御供や人体損壊含む)の習慣、始まる。 1年 - Wikipedia い(1)けにえの習慣マヤではじまった。 西暦2年 仏教が中国に伝わる。 2年 - Wikipedia 中国に(2) 仏教きたよ。ありがたや。 西暦3年 野見宿禰が埴輪を考案。 3年 - Wikip…

かけ算九九の表をグラフにしてみました

library(dplyr) library(cowplot) kuku <-data_frame(left=rep(1:9,each=9),right=rep(1:9,9)) %>% mutate(ans=left*right) ggplot(kuku,aes(x=factor(left),y=factor(right)))+ geom_tile(aes(fill=ans))+ geom_text(aes(label=ans))+ scale_fill_continuou…

エイプリルフール

今日がエイプリルフールだからってぼくの言葉にはひとかけらの嘘もない。無駄なく、よどみなく、取り替えがきかない、まるで俳句のように、安心しろすべてがうまくいく。

pymc3 で状態空間(ローカルレベル)モデル

Python ハローワールドです。Rodeo(Yhat End-to-End Data Science Platform: Rodeo)っていう RStudio みたいなやつをいれました。なかなか快適です。 import numpy as np import pymc3 as pm import matplotlib.pyplot as plt import statistics as stat i…

R によるすごくかんたんなパーティクルフィルタの実装例

これであってるのかあんまり自信ない。主に http://daweb.ism.ac.jp/koza/koza2008/PF_Nakano20081030.pdf を参考にしました。パーティクルフィルタのアルゴリズムは、 一期先の予測を乱数でばらまく 得られた乱数に尤度の重みをつけてまき直す だと理解しま…

Rcpp で独立メトロポリスヘイスティングス

R

独立メトロポリス・ヘイスティングス法を用いたベイズ推測の簡単な例題 - 廿TT でやったのと同じことを Rcpp で書いてみた。 ハローワールド。C++ のコードはこう。 #include <Rcpp.h> using namespace Rcpp; double lik(double lambda, NumericVector dat) { double</rcpp.h>…

[searchConsoleR]テキストマイニングことはじめ:検索キーワードの視覚化(ワードクラウド、ワードカウント、共起ネットワーク)

R SEO

はじめに フリーソフトだけでテキストマイニングしたい。R と MeCab を使います。 MeCab: Yet Another Part-of-Speech and Morphological Analyzer RMeCab - RとLinuxと... R のパッケージもいっぱい使います。 library("googleAuthR") library("searchConso…

未知の変化点があるモデルでは AIC が使えない

モデル 時系列データ () があるとします. このデータが, 変化点()以前では平均 , 標準偏差 1 の正規分布に従い, 変化点から後には平均 , 標準偏差 1 の正規分布に従うと考えます. 標準偏差は既知とします.ここで は標準正規分布に従う確率変数です.変化点 …

ゼロ過剰負の二項分布によるセッションの間隔のモデル(Google アナリティクス)

Google アナリティクスではセッションの間隔(daysSinceLastSession)という指標を見ることができる.これはサイトに訪れたユーザーの直前のセッションが何日前だったのかを示すものだ.ご覧の通り非常にゼロの多いデータで, しかもロングテールなので, これを…

ggplot2 で欠けた円グラフ

久しぶりに円グラフを描きたくなった。100% に満たない量を表す欠けた円グラフ。 library(cowplot) dat<-data.frame(group=LETTERS[1:2],y=c(0.7,0.4)) ggplot(dat)+ geom_bar(aes(x=group,y=y*100),width=100,stat = "identity")+ ylim(c(0,100))+ facet_wr…

ガンマ再生過程に基づくカウントデータの分布

関心のある事象(例えば機械の故障, タクシーの到着など)が繰り返し生起し, それぞれのイベントの生起間隔が独立に同一のガンマ分布に従う場合を考える.イベントの生起間隔を確率変数 で表す. またイベントの発生時刻は, で表す. いま, 開区間 で起こったイ…

n 人をふたつにわける組み合わせを列挙する

R

ふたつの組 A と B があり, を i 番目の人を組 A に入れるとき 1, 入れないとき 0 を取る変数とする. 下のような樹形図で考えると, 2 値の変数が n 個あるので, 組み合わせの数の総数は . (Binary Tree clip art Free Vector / 4Vector より)なので最初に …

Google アナリティクスのインタレストカテゴリを平行座標プロットで再クラスタリング

インタレストカテゴリとは Google ではオンラインでの活動や購買行動からユーザーの興味・関心を推測して、ユーザーを分類しています。この分類は「インタレストカテゴリ」と呼ばれています。Google アナリティクスでは、インタレストカテゴリには「アフィニ…

非定常ポアソン過程でアフィリエイトのコンバージョンを予測(グループドデータ版)

abrahamcow.hatenablog.com の続きです。めったに起きないコンバージョンの成長の非定常ポアソン過程によるモデル - 廿TT では一日に複数のコンバージョンが発生しても、それは一回とカウントしていました。これはもったいない。こういうのはグループドデー…

日付データを年ごとや月ごとや週ごとで集計

R

R です。lubridate パッケージを使うとかんたん。[1301 東証1部] 極洋 日足 時系列データ CSVダウンロード のデータを使います。 > library(lubridate) > library(dplyr) > stocks <- read.csv("~/Downloads/stocks_1301-T.csv",fileEncoding = "cp932") > h…

AWK でかけ算九九の表を作る

AWK

$ seq 9 | awk '{for(i=1;i<10;i++)$i=$1*i;print}' 1 2 3 4 5 6 7 8 9 2 4 6 8 10 12 14 16 18 3 6 9 12 15 18 21 24 27 4 8 12 16 20 24 28 32 36 5 10 15 20 25 30 35 40 45 6 12 18 24 30 36 42 48 54 7 14 21 28 35 42 49 56 63 8 16 24 32 40 48 56 64…

状態空間非定常ポアソン(NHPP using Stan)

ポアソン過程は再発事象のモデルとしてよく使われる。ポアソン過程ではイベントが観測された時刻を () とすると、イベントの生起間隔 は独立にパラメータ λ の指数分布に従う。ポアソン過程の拡張としてパラメータλ が時間に依存して変化する非定常ポアソン…

(R+Google アナリティクス)エラーバーで信頼下限をプロット

場面設定 当サイトは女性の訪問者が少ないので、女性の訪問者を増やしたいと思っている。サイト作りの参考にしようと女性の新規訪問の割合が多いランディングページをリストアップしたい。そこで新規セッション率で降順にソートをかけると、セッション数10〜…

ラザニアプロット(fields パッケージの image.plot にちょっと一工夫)

このエントリは計算機統計学会第30回シンポジウムにおける兼田麻里奈、坂本亘両氏のご発表「ラザニアプロットを用いた経時データの視覚化」にインスパイアされたものです。時系列データを表現するのによく用いられるのは折れ線グラフですが、複数の系列を同…

rstan で混合二項分布のパラメータ推定

ordered 型は、「小さい順」という制約です。StanとRでベイズ統計モデリングで解説されている「ラベルスイッチング」を回避するためにこれを使ってます。StanとRでベイズ統計モデリング (Wonderful R)作者: 松浦健太郎,石田基広出版社/メーカー: 共立出版発…

Missing Not At Random(MNAR):R と Stan で欠測が欠測データに依存する場合のパラメータ推定

測定機器かなにかの都合上、観察対象の値が小さくなると欠損が出やすくなる状況を考えます。R で以下のようにしてシミュレーションデータを生成しました。 set.seed(1) N <-200 X<-rnorm(N,2) X2 <-ifelse(runif(N)

R {deSolve} で差分方程式

R の deSolve パッケージで差分方程式を計算するには ode のメソッド "iteration" を使う。関数 func は変化率ではなく状態変数の新しい値を返すように書く。 library(deSolve) Ti <- 100 a <- 1.2 disc_logis <-function(Time,x,a){ x2 =a*x*(1-x) list(x2)…

Stan の integrate_ode_rk45 を使ってバスモデルのパラメータ推定

以下の微分方程式で記述されるモデルのパラメータを推定します。ほんとうは閉じた形で解が求まるのですが今回は Stan の integrate_ode_rk45 を使って数値的に解を求めます。カラーテレビの普及率のデータ(第1章第2節3 1 情報通信機器の世帯普及率 : 平成1…

AWK で Reservoir Sampling; テキストからランダムに少数の行を抽出

AWK R

R による溜池サンプリング(Reservoir Sampling)の実験 - 廿TT を踏まえて, AWK でテキストファイルからランダムに1000行非復元抽出するコードを書きました.テスト用のデータをRで生成します. set.seed(1) rmixnorm3 <- function(n) { n1 <- round(n*0.5) n…

R による溜池サンプリング(Reservoir Sampling)の実験

R

Data Stream Management という本に出てくる Reservoir Sampling(溜池サンプリング)という手法をシミュレーションしてみたい.これはサイズ N の母集団(N は未知でもよい)からサイズ n()のランダム標本を非復元抽出で取ってくるアルゴリズムで, 大きす…

R {Nippon} パッケージのコロプレス図(塗り分け地図)に凡例をつける

全国最低賃金 地域別最低賃金の全国一覧 |厚生労働省 の表をプロットしてみる.カラーパレットには RColorBrewer パッケージを使います. library(rvest) library(dplyr) url1 <-"http://www.mhlw.go.jp/stf/seisakunitsuite/bunya/koyou_roudou/roudoukijun/…

負の二項分布を用いたリピート回数のモデル(Google アナリティクス)

序 負の二項分布は「試行が n 回成功するまでに失敗した回数 x の分布」として知られている( 負の二項分布 | r回の成功を得るのに必要な試行回数 )一方で, ポアソン分布のパラメータ λ がガンマ分布するような分布でもある( 可視化で理解する「負の二項分…

R {arules} によるアソシエーション分析の結果をどうやったら見やすく表示できるか試行錯誤中

library(arules) library(cowplot) data("Adult") rules <- apriori(Adult, parameter = list(supp = 0.5,conf = 0.9,target = "rules",maxlen=2)) rules_lhs <-as(lhs(rules),"list") #条件部 rules_rhs <-as(rhs(rules),"list") #結論部 rules_lhs <- sapp…

R {mixtools} で混合多項分布のパラメータ推定をやってみた

R

場面設定 全く同意できない 同意できない 同意できる 非常に同意できる のような四択形式のアンケート20問を100人に対して行ったとする. 問題ごとに質問項目への回答数を集計してプロットしたのが以下の図である. アンケートの反応パターンから, 質問項目を…

階層ベイズでもサンプルサイズを増やしたらベイズ信頼区間の幅は細くなってくれるのか

系列長 の時系列データがあるとして, これに以下のようなモデルを当てはめます. 階層事前分布として , には平均 10 の指数分布, には幅の広い一様分布を仮定します. いま の95%ベイズ信頼区間を求めたいとします.ベイズ信頼区間の幅はサンプルサイズ(系列長…

EMアルゴリズムの練習:右打ち切りデータからの指数分布のパラメータ推定

アルゴリズム 独立に平均 の指数分布に従う大きさ の標本 があり, 内 は完全な観測が得られ, は右打ち切りされているとする. 完全な観測が得られたとして, このときの対数尤度は, である. ただし, は観測されない打ち切られた線分の長さを拡大して表したもの…

左切断指数分布の期待値

が平均 の指数分布に従うとする. という条件付きの密度関数は以下で与えられる. 期待値を取ると, 弱点克服大学生の確率・統計作者: 藤田岳彦出版社/メーカー: 東京図書発売日: 2010/04/09メディア: 単行本購入: 6人 クリック: 9回この商品を含むブログを見る