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

廿TT

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

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

ロトカ・ヴォルテラの方程式については、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.1.0)

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 で前回訪問…

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

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

[searchConsoleR]要対策キーワードをあぶり出す回帰分析

R SEO

辻正浩氏によると、 ある程度以上の表示回数があり、かつ、ある程度以上の掲載順位があるが、クリック数がある程度以上になっていないキーワードを抽出するだけで、タイトル要素に何か問題があるのかもしれないとか、検索されている意図と実際のサイトが合っ…

(R + Google アナリティクス)ランディングページの集客力推移のヒートマップ

多くの訪問を集める主要ページにはこういう施策、テールの部分を支える裾野ページにはこういう施策、といった具合に、ランディングページを層別にして対策を立てる場合を考えます。このとき、訪問(セッション)数全体の推移を見ただけでは、どの層がどう変…

twitteR でかんたんレピュテーション・マネジメント(評判分析)

レピュテーション・マネジメントとは twitteR が使えるようになった(twitteR を使えるようにするためのメモ(version 1.1.9) - 廿TT)のでかんたんな集計を行ってみます。レピュテーション・マネジメントとは企業や組織がみずからの評判を守るために行う活…

twitteR を使えるようにするためのメモ(version 1.1.9)

R から twitter のデータを引っ張ってこれるパッケージ twitteR、だいぶかんたんになってました。まず最初にじぶん用のツイッターアプリをつくる必要があります。https://dev.twitter.com/ にツイッターアカウントでログイン。下のほうにちっちゃく出てる Ma…

ggplot2 で Google アナリティクスデータの人口ピラミッド

下表のデータを人口ピラミッドの形で描いてみる。(もうちょっときれいに描ける人いそう。ご意見求む。) usergender userage bracket users female 18-24 628 female 25-34 1001 female 35-44 416 female 45-54 174 female 55-64 85 female 65+ 40 male 18-…

(R + Google アナリティクス)2014年人気だった記事の集客力はどのくらい目減りしたか

前置き コンテンツをつくるのはお金や手間がかかるので、ひとつひとつの記事が集客力を維持してくれるとうれしい。そこで2014年のセッション数トップ10の記事が、2015年にはどのくらいセッション数があったか集計してみる。集計には RGA パッケージを使う。R…

変化点のあるポアソン過程のパラメータの最尤推定

尤度関数 変化点のあるポアソン分布のパラメータの最尤推定 - 廿TT では、生起したイベントの個数に着目しましたが、生起の間隔に着目してモデル化することもできます。 とすると、変化点のない強度(intensity)λ の定常ポアソン過程では、点と点の間隔は指…

RGoogleAnalytics の使い方(Version 0.1.1)

R から API で Google アナリティクスのデータを引っ張ってくるパッケージ、RGoogleAnalytics を使えるようにするまでのメモです。パッケージのインストールは install.packages("RGoogleAnalytics") でいけます。まずは https://console.developers.google.…

Nerlove-Arrow の広告−販売モデル

モデル Nerlove と Arrow(読み方はネルロフとアロウでいいのかな?)が提案したらしい広告に対する市場の販売のモデルというのがあり、これは以下の通り教科書に出てくる「一階線形常微分方程式」そのものの形をしている。 A(t) が時刻 t での売上 q(t) が…

変化点のあるポアソン分布のパラメータの最尤推定

モデル 変化点のあるポアソン分布についてはいろんな研究がなされていますが、一番単純と思われる手法を試します。解析対称はイギリスの炭鉱事故の発生件数のデータです。3. Tutorial — PyMC 2.3.6 documentation から取得しました。 #R のコード dat <- c(4…

めったに起きないコンバージョンの成長の非定常ポアソン過程によるモデル

弱小アフィブログの成長 これは当ブログのアマゾンアフィリエイトレポートの折れ線グラフです。横軸が日付, 縦軸がその日の注文数合計です。注文数合計は 0 が多く, 折れ線グラフは下に潰れて, 注文があった日だけスパイクのようにとげとげつきだしています…

観測期間にギャップのあるワイブル過程のパラメータの最尤推定

最尤推定量 ワイブル過程のパラメータの最尤推定 - 廿TT の続きです.Crow (1988) に出てくる例題をやります. を区間 でのイベントの発生時刻, を区間 でのイベントの発生時刻として,非定常ポアソン過程を区間 と ()で観測した場合の尤度関数は以下のように…

ワイブル過程のパラメータの最尤推定

ワイブル過程とは ワイブル過程とは強度(intensity)関数に を仮定した非定常ポアソンのことで, おもに信頼度成長のモデルに使われる.この強度関数はワイブル分布のハザード関数と一致する.ワイブル分布とは別ものなので紛らわしい.power-law ポアソン過程…

RStanで離散パラメータを含むモデルの推定(disaster model)

前置き Stan のマニュアル11章の例題です。 Stan は離散パラメータをサポートしていないので、離散パラメータを含むモデルの推定では周辺化して離散パラメータを消去してやる必要があります。その練習です。解析対称はイギリスの炭鉱事故の発生件数のデータ…

最尤法を理解するためのお絵かき

R

黒い点が観測値 線分で表した高さの積が尤度 赤と青では青のほうが尤度が大きい パラメータを動かして最も尤度が大きい分布を見つけるのが最尤法 最尤法の説明はむずかしくて、下手すると半端に哲学っぽくなっちゃうんだけど、ぼくの場合はこういう図を描い…

超球の体積をモンテカルロシミュレーションで求める

R

モンテカルロ法の例題として、円周率πの値の近似はよく見る。正方形の領域に一様に点を打ち、内接する円の中に入った点の個数を全部の点の個数で割ってやると、正方形の面積と円の面積の比が出てくる。正方形の面積と円の面積の比に正方形の面積を掛けてやる…

RStan でベイズ逆問題もどき

, という系列を考える。A は定数。(この系列になんらかの解釈があればおもしろかったんだけど思いつかなかった。)いま、観測されるのは のみである。 が既知ならば A を推定するのはかんたん。A が既知ならば を推定するのはかんたん。しかし、いま観測さ…

離散時間データからのワイブル分布のパラメータの最尤推定

尤度 イベント発生が区間 に起こったことがわかっており、イベントが発生した時間そのものはわからない状況を考える。このようなデータを区間打ち切り(interval censored)データと呼ぶ。尤度は、である。 シミュレーション R の survival パッケージでは T…

離散時間データからの指数分布のパラメータの最尤推定

尤度関数と最尤推定量 イベント発生が区間 に起こったことがわかっており、イベントが発生した時間そのものはわからない状況を考える。このようなデータを区間打ち切り(interval censored)データと呼ぶ。尤度は、である。以降、指数分布を仮定して考える。…

[R]度数分布表を作って折れ線グラフで比較する

set.seed(1234) X <- rnorm(100,1,1) Y <- rnorm(100,3,1) Z <- rnorm(100,2,0.5) breaks <-seq(-5,7,by=1) plot(table(cut(Z,breaks)),type="b",col="royalblue",ylab="frequency") lines(table(cut(Y,breaks)),type="b",col="forestgreen") lines(table(cu…

独立メトロポリス・ヘイスティングス法を用いたベイズ推測の簡単な例題

R

なぜ MCMC が必要か ベイズ統計学では母数 θ を確率変数と考える。ベイズ統計のための準備, ベイズの定理, 事前分布と事後分布 - 廿TTベイズ推測はデータ が得られたもとでの θ の確率分布 事後分布を通じて行われる。母数 θ の点推定には事後期待値(expect…

一階線形常微分方程式を使った流行語のモデル

モデル化 流行語は話題にする人が増えれば増えるほどさらに流行しやすくなり、ある時期を超えると忘れられていくと考えられる。このダイナミクスをなるべく簡単な微分方程式で記述したいのでこのような仮説を置いた。 流行語の拡散スピードは話題にされる量…

タンクモデルのシミュレーション

タンクモデルとは 気象庁|土壌雨量指数 気象庁|流域雨量指数 を参照。 (上記ページから引用)シンプルといえばシンプルだけどこんなモデルどうやってパラメータ推定するんだろう。 都市域:1段タンクモデル #R のコード library(deSolve) library(tidyr) …

RStan でポアソン過程の停止時刻を推定する

問題設定 n = 20 人の患者さんが繰り返し訪問する病院の窓口を T = 100 時間観察した。患者 i はそれぞれある時期 に達すると病院への訪問をやめる(理由は不明:病気が治ったか、遠くに引っ越したか、通院に飽きたか、死んだか)。この訪問がレート λ (全…