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

廿TT

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

[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回この商品を含むブログを見る

(ggplot2)順位の入れ替わりプロット

なんと呼ぶのかわからないグラフ library(cowplot) set.seed(100) old <- data.frame(Group = "old", Rank = 1:5, Text = c("一郎","次郎","三郎","四郎","五郎")) new <- data.frame(Group = "new", Rank = sample(1:5), Text = c("一郎","次郎","三郎","四…

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

下記のブログで公開されている拡張カルマンフィルタのコードでもうちょっと遊んでみる. 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 のコードは以下のようになります。 # 後…

ツーペン(Toepen)というゲームアプリを買った

おもしろいのでアフィリエイト広告経由でみんなも買ってくれ。ToepenGraafICTゲーム¥150ToepenHD - het leukste kaartspel Toepen op je iPad!GraafICTゲーム¥250ゲームの遊びかたは以下に書きます。 大まかな流れ speel でゲームスタート。ツーペン(Toepen…

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

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

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

R

最尤推定と条件付き最尤推定 ゼロ切断ポアソン分布のパラメータの最尤推定 - 廿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, the free encyclopedia 最尤推定量 サンプルサイズを n とすると対…

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

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

AWK で gather と spread

AWK R

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

AWK で新しいことわざを作る

AWK

動機 カットアップで小説を書いたウィリアム・バロウズという作家をご存知だろうか。カットアップとは印刷物の文章をハサミで切って並べ替え、おもしろそうなフレーズができたらそれをそのまま使うという手法だ。実にお手軽。誰にでもまねできそうだ。しかし…

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…