廿TT

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

R

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アルゴリズムの練習:右打ち切りデータからの指数分布のパラメータ推定

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

拡張カルマンフィルタ:ロジスティック方程式

下記のブログで公開されている拡張カルマンフィルタのコードでもうちょっと遊んでみる. R code for implementing the extended Kalman filter – R code, simulations, and modeling R code for estimating the parameters of an extended Kalman filter mode…

拡張カルマンフィルタによるロトカ・ヴォルテラ方程式へのデータ同化

※ぼくは「データ同化」のなんたるかについてはまるでわかっていない.下記のブログで解説されている, 拡張カルマンフィルタを試してみる. R code for implementing the extended Kalman filter – R code, simulations, and modeling R code for estimating th…

R によるシミュレーテッド・アニーリング法の練習

アルゴリズム 関数 を最大にする を探す問題を最適化問題という. シミュレーテッド・アニーリング法は以下のようなアルゴリズムで最適化問題を解く. 反復 t で, をシミュレートする. 確率 で を受理する. さもなくば とする. ここで . また密度 g は原点 0 …

一階線形常微分方程式の解のピークの位置

計算 一階線形常微分方程式とは以下の形で表される微分方程式のことである. 解は, が最大値をとる点を求めるには, と置いて解けばよい. 部分積分の公式を逆に使って, この方程式を解けば, が最大値をとる点が求まる. ここで,より, つまり が最大値をとる点を…

一次元の拡散方程式をノイマン境界条件の下で解く(差分法、陰解法)

陰解法 以下の偏微分方程式が拡散方程式である。ここで x は空間方向の位置を表し、t は時間を表す変数である。 が求めたい未知関数で、初期条件 は与えられた関数である。D は正の定数である。条件 、 のように端点における微分の値を固定する条…

打者の調子の波のモデル化 4

データ MLB - スポーツナビ にプロ野球選手ボガーツの打席結果が掲載されています。例として下表に一部を抜き出します。 日付 打数 安打 6/1 5 1 6/2 5 1 6/3 5 2 6/4 3 0 6/5 4 3 6/6 4 0 このデータをとってくる R のコードは以下のようになります。 # 後…

ゼロ切断ポアソン分布によるセッション数のモデル化とデータ拡大による潜在利用者数の推定

序 ゼロ切断ポアソン分布によるセッション数のモデル - 廿TT ではウェブサイトのセッション数をゼロ切断ポアソン分布としてモデル化し、ウェブサイトにアクセスしなかった人も含めた潜在的な利用者数 N を推定しようとした。しかし日によって N の推定値が大…

EM アルゴリズムによるゼロ切断ポアソン分布のパラメータ推定

最尤推定と条件付き最尤推定 ゼロ切断ポアソン分布のパラメータの最尤推定 - 廿TT の続きです。上記のエントリでは、0 でない観測の数 n を所与として、n 個の観測が得られるために必要な試行回数 N は観測されないとしていました。しかし現実的には n が固…

阿部誠(2008)RF 指標から生存確率を求める(rstan 版)

モデルと尤度関数 阿部(2004)RF指標から生存確率を求める - 廿TT の続きです。阿部(2004)によれば, 以下の 2 つの仮定を置くことで, RF(リセンシーとフリクエンシー)指標のみから, 顧客の生存期間を求めることができたのでした.フリクエンシーを x, リ…

ゼロ切断ポアソン分布によるセッション数のモデル

あるウェブサイトの利用者が全部で N 人いるとして、N 人がある一定時間内に 回そのウェブサイトを利用するとします。一定時間内に一回以上サイトを訪れた利用者の数を n とすると、 n に相当するのがユーザー数です。一定時間内にサイトを訪れるユーザーの…

ゼロ切断ポアソン分布のパラメータの最尤推定

確率質量関数 f をポアソン分布の確率質量関数とすると、ゼロ切断ポアソン分布(zero-truncated Poisson)の確率質量関数は、期待値は、Zero-truncated Poisson distribution - Wikipedia 最尤推定量 サンプルサイズを n とすると対数尤度関数は、これを最大…

時系列で見て安定している KPI の探し方

要約 見よう見まねでぼくのブログの KPI を策定してみた。 折れ線グラフでスパイクをチェック 変動係数で散らばりをチェック 箱ひげ図で曜日効果をチェック 良き KPI とは KPI は時系列で見て安定していることが求められる。tokoroten 氏は KPI を乱高下させ…

AWK で gather と spread

AWK R

はじめに dplyr ユーザーのための AWK 入門 - 廿TT に引き続き、相模原市オープンデータライブラリー | 相模原市 で公開されている駅別乗降人員の推移データを使用して AWK で簡単なデータ整形を行ってみます。 列持ちのデータを行持ちに変える wide フォー…

dplyr ユーザーのための AWK 入門

R AWK

はじめに dplyrを使いこなす!基礎編 - Qiita を参考に、相模原市オープンデータライブラリー | 相模原市 で公開されている駅別乗降人員の推移データを使用して dplyr と対比させながら AWK で簡単な集計を行ってみます。 行の絞り込み dplyr でいう filter …

R: 複数回答のアンケートのグラフの例

経緯 「新知事に優先して取り組んでほしい政策は「教育・子育て」が41・9%で最多」都知事に改憲阻止を求めている有権者などいないようである (゚⊿゚) 【東京都知事選 序盤情勢】https://t.co/TheBuRLRPV pic.twitter.com/NOWyQFxQSb— ano_ano (@ano_ano_an…

時系列データで相関を出してはいけないのか(失業と自殺は関係あるのか2)

経緯 じ、時系列データに対して単純な相関を算出している。。。 https://t.co/3yUB5ZEhRo— 統計たん@Rアイドル (@stattan) 2016年7月15日男性に関して言えば、失業率と自殺率は強い相関を持つことが舞田敏彦らによって指摘されている。相関係数は0.7224。デ…

Stan コードをベクトル化したら速くなるか:AR(1) モデルを例に

ちょっと速くなる。こういうモデルのパラメータを推定しました。 は平均 0、分散 の正規分布、 の事前分布は一様分布です。ベクトル化前のコードはこちら。 // AR1_1.stan data { int<lower=1> T; real y[T]; } parameters { real beta0; real beta1; real<lower=0> sigma; } m</lower=0></lower=1>…

[KFAS]0-1データの状態空間モデル(打者の調子の波のモデル化 3)

とりあえずプロット Albert (2008) 打者の調子の波のモデル化 - 廿TT 打者の調子の波のモデル化(幾何分布編) - 廿TT と同じく、カルロス・ギーエンの打撃成績のデータを使います。 y <- c(0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,…

ロトカ・ヴォルテラの方程式のパラメータ推定

ロトカ・ヴォルテラの方程式については、abrahamcow.hatenablog.comを見てください。 hare–lynx データ ここに一組のデータ・セットがあります。 Hare=c(30, 47.2, 70.2, 77.4, 36.3, 20.6, 18.1, 21.4, 22, 25.4, 27.1, 40.3, 57, 76.6, 52.3, 19.5, 11.2, …

ggplot2 で接線場や方向場を描く

接線場 微分方程式の解の振る舞いをみるために接線場を描いてみる。以下で与えられるロジスティック方程式を考える。接線場は t-x 平面の各点に傾き の小さな線分を描いたものである。 logis <- function(t,x,a=1){ a*x*(1-x) } #点(t1,x1)を通り, 傾きが sl…

自殺者数の推移(原因・年代別)

自殺白書に関する報道は東京新聞と産経新聞で見出しの付け方が対照的だった。 東京新聞:若者と高齢者の自殺深刻 政府、16年版白書:社会(TOKYO Web) 自殺白書 「経済」理由の自殺半減 厚労省「法改正や法律相談の充実が奏功」 (産経新聞) - Yahoo!ニュー…

Alluvial diagram で Google アナリティクスのデータをプロットして遊ぼう

Alluvial diagram(アルビアルダイアグラム)をかんたんにプロットできる ggalluvial パッケージが、Rで解析:Alluvial diagramsをプロットしませんか「ggalluvial」パッケージ で紹介されていた。これを使ってみようと思う。 入口ページと出口ページ まずは…

googleAnalyticsR の使い方(Version:0.4.0)

最終更新日:2017年6月14日 V4 API の概要 Google Analytics Core Reporting API V4 では以下のような機能が追加されました。 データの範囲を 2 つ指定できるようになった。 指標を四則演算して新たな指標を作り出すことができるようになった。 コホートレポ…

[rstan]横浜駅SFの移入と離脱のモデル

横浜駅SF (カドカワBOOKS)作者: 柞刈湯葉,田中達之出版社/メーカー: KADOKAWA発売日: 2016/12/24メディア: 単行本この商品を含むブログ (4件) を見る 経緯 この減衰曲線はどういう関数で説明するのが適切か考えている(指数関数ではない)https://t.co/PYdRO…

気温・湿度と熱中症の危険度(暑さ指数)の表

厚生労働省:職場における熱中症の予防について より暑さ指数(WBGT)の表をプロットしてみた。 WBGT値 警戒水準 25℃未満 注意 25℃~28℃ 警戒 28℃~31℃ 厳重警戒 31℃以上 危 険 そろそろ熱中症に気をつけましょう。 library(pdftools) library(pipeR) librar…

岩波データサイエンス Vol.1 の年輪の例題を dlm でやる

DLM ノイズが正規分布し、変数間の関係が線形の状態空間モデルは動的線形モデル(Dynamic Linear Model; DLM)と呼ばれる。動的線形モデルは以下の式で表現できる。 推定問題 現時刻を k とするとき、推定問題はつぎの三つに分類できる。 予測(prediction)…

[dlm]状態空間モデルでトレンドと広告の効果を分離して推定する

はじめに Stanで統計モデリングを学ぶ(7): 時系列の「トレンド」を目視ではなくきちんと統計的に推定する - 東京で働くデータサイエンティストのブログ をみてください。上記事では Stan で状態空間モデルを推定しているので、ここでは R の dlm パッケージ…

左切断右打ち切りデータからの最尤推定

切断と打ち切り 観測対称がとる値がある範囲を超えたとき、まったく観測されない場合を切断(truncated)と呼ぶ。それに対して、観測されない対象の個数はわかっている場合を打ち切り(censored)と呼ぶ。観測範囲がある値以上に限定されている場合、左切断…

Google アナリティクスデータから指標を偏差値化して記事を評価

偏差値好きな人多い印象あるけどぼくはあんまり好きじゃなかった。でもなんとなく Google アナリティクスデータで偏差値出してみたら、これはこれでけっこういいかも、と思った。次元(単位)のない量にして、複数の指標どうしを概観的に比較できる。セッシ…

[searchConsoleR]検索クエリとページタイトルのマッチングをはかる

SEO R

検索結果画面に表示されるページタイトルが検索クエリとマッチしていない場合、そのページはあまりクリックしてもらえない。では検索クエリとページタイトルのマッチさせるにはどうするか。これはむずかしい。考えられるもっともシンプルな方法は、検索クエ…

大森・宇津公式を強度関数とした非定常ポアソン過程の重ね合わせによる熊本地震の余震の分析

スクレイピング 使用するツールは R です。データは 気象庁|地震情報 から取得しました。rvest パッケージを使うとかんたんです。また XML パッケージの readHTMLTable 関数を使う方法もあります。R から HTML の表を読み込む - 廿TT速報は一定時間たつと流…

MASS::fitdistr で打ち切りデータからの最尤推定

R には fitdistr という最尤推定をしてくれる関数がありまして、自分で密度を定義してやれば打ち切りデータからでもパラメータを推定できる。 mydens <-function(x,shape,scale){ ifelse(d==1, pweibull(x,shape, scale,lower.tail = FALSE), dweibull(x,sha…

エクセルで無相関検定:失業率と野菜摂取量の相関

町山智浩がアメリカでは貧しい人たちが野菜を食べれなくて困っているというような話をしていた。町山智浩 映画『Fed Up』が描くアメリカの飢餓・肥満問題を語る日本でもそういう現象があるだろうと思い、野菜摂取量と失業率の相関を調べてみた。使用したデー…

打者の調子の波のモデル化(幾何分布編)

仮説検定 Albert (2008) "Streaky Hitting in Baseball" ではベータ二項分布を用いて野球選手の調子の波を評価した。Albert (2008) 打者の調子の波のモデル化 - 廿TT下記はカルロス・ギーエンという選手の2005年の打撃成績のデータで、ヒットを 1、アウトを …

optim による glm:最尤推定 + 信頼区間, Wald 検定

あらまし 自分で尤度を書いてみて, R の glm 関数がやってることを再現する. 尤度さえ書ければパラメータの点推定, 区間推定ができるし検定もできる. それができるようになれば, パッケージなどが用意されていない新しいモデルでも計算できるようになる(と…

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

モチベーション 野球にはまったく興味ないんだけど、Albert (2008) "Streaky Hitting in Baseball" を読んでた。http://citeseerx.ist.psu.edu/viewdoc/download?rep=rep1&type=pdf&doi=10.1.1.150.5808スポーツデータ解析に関する論文を探すには Journal of…

新規率と直帰率の関係からランディングページの改善点を探る

はじめに [対談]SEO 辻正浩氏×アクセス解析 小川卓氏:検索キーワードが見えない時代のSEOとは? | Web担当者Forum で紹介されていた 検索結果における各ページのクリック数と掲載順位の関係をグラフにすると1つの曲線に集約される。曲線から大きく外れた…

阿部(2004)RF指標から生存確率を求める

R

はじめに 阿部誠(2004)の手法をシミュレーションで再現してみます. http://merc.e.u-tokyo.ac.jp/mmrc/dp/pdf/MMRC16_2004.pdf顧客関係管理(CRM)の現場では、RFM 分析という手法が広く使われているそうです. 顧客の購買行動の詳細を蓄積していない企業で…

多項分布のパラメータ推定:モンテカルロ EM アルゴリズムの例題

R

序 『R によるモンテカルロ法入門』p.177 の例題をやります. 『R によるモンテカルロ法入門』は書き方に省略が多くて意味がわからないところがけっこうあるので『計算統計学の方法』pp.82-pp.88 もあわせて参照しています.Rによるモンテカルロ法入門作者: C.…

犯罪の検挙率のグラフ(モザイクプロット)

犯罪統計資料(統計表一覧 政府統計の総合窓口 GL08020103)によると、平成28年1~2月の検挙率は35.2%でした。過半数の犯罪は検挙されていません。しかしこの数字はシンプソンのパラドックスめいたところがあって、犯罪のうち、凶悪犯、粗暴犯、知能犯、風…

貸借対照表のグラフ(コロプラ)

特に理由はないけれど、コロプラの貸借対照表を図示してみた。貸借対照表 | 業績・財務 | IR情報 | 株式会社コロプラ library(readxl) dat_row <-read_excel("~/Downloads/colopl_financialdata_201512.xlsx",9,skip=1) tmp <-dat_row[,-c(1,3:4)] colnames(…

Google アナリティクスデータの変化を検知するための CUSUM 管理図

問題 例えばソーシャルメディアで炎上してたりしたらサイトの訪問(セッション)数は伸びるだろう。セッション数を定点観測していればなんらかの異常があったときにすぐに対処できる……はずだ。Googleアナリティクスオタクの私が、毎日見ているたった1つのレ…

(R + Google アナリティクス)継続率を計算してみる

継続率にはいろんな定義がありえる。ここで計算するのは「ある日に初回訪問したユーザーのうち、x日以内に再訪問したユーザーの割合」のこと。Google アナリティクスのディメンション、sessionCount で訪問の回数がわかる。daysSinceLastSession で前回訪問…

状態空間モデルで自然検索トラフィックの成長を予測する

場面設定 コンテンツを増やせばそれだけ自然検索にヒットするページが増え、ウェブサイトのトラフィックは増加します。向こう一年間これだけ記事を書くぞ、というのが決まっていたとして、その計画から自然検索経由の訪問(セッション)数を予測できるでしょ…