廿TT

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

R

区間打ち切りデータからの新型コロナウイルスの潜伏期間の推定(rstan版)

背景の整理 COVID-19 の潜伏期間をrstanで推定する - 驚異のアニヲタ社会復帰の予備 で紹介されている論文 Incubation period of 2019 novel coronavirus (2019-nCoV) infections among travellers from Wuhan, China, 20-28 January 2020. - PubMed - NCBI …

Turnbull のアルゴリズム(R の survival パッケージ)を用いた新型コロナウイルスの潜伏期間の推定

COVID-19 の潜伏期間をrstanで推定する - 驚異のアニヲタ社会復帰の予備 の紹介で知った。 https://docs.google.com/spreadsheets/d/1jS24DjSPVWa4iuxuD4OAXrE3QeI8c9BC1hSlqr-NMiU/edit#gid=1449891965 に,感染者が新型コロナウイルスのある地域への暴露を…

状態空間トピックモデルで詩の流れを可視化

中原中也のサーカスという詩をたぶんあなたはすでにご存知だろう. サーカス: 中原中也・全詩アーカイブ 「幾時代かがありまして」の反復ではじまり「ゆあーん ゆよーん ゆやゆよん」の反復で終わる. この流れを可視化したい. モデル をワン・ホット・エンコ…

メタゲノムデータとメタボロームデータで相関を出してはいけないか

R

ちょっと質問があったので R でシミュレーションしてみます。 とはいえこの説明でふつうのお医者さんが理解できるかは疑問。 背景の整理 メタゲノムデータとは サンプルに含まれるDNAをまとめて回収してきて、そのDNAの配列を読み取り、データベースに照会し…

ggplot2 などのモダンなパッケージを使わずにトピックモデルのパラメータを可視化したい

腸内細菌のデータを使います。 説明は気が向いたら書きます。 library(curatedMetagenomicData) plot.minibarTable <- function(tab,layoutMat,leftmargin=15, col="black"){ oldpar <- graphics::par(no.readonly = TRUE) N <- length(tab) ran <-c(0,max(s…

geom_pointrangeみたいなやつを並べて書きたい

set.seed(1) meanse <- lapply(1:7, function(i)cbind(mean=rnorm(10,0,2),se=rexp(10))) names <- sapply(1:10,function(i)paste0(sample(LETTERS,10),collapse = "")) tab <- list(meanse=meanse,names=names) class(tab) <- "seTable" plot.seTable <- fu…

棒グラフをとにかくいっぱい並べて書きたい

set.seed(2) y <- lapply(1:7, rexp, n=26) names(y) <- sapply(1:7, function(x)paste0(sample(LETTERS,10),collapse = "")) tab <- list(y=y,x=LETTERS) class(tab) <- "minibarTable" plot.minibarTable <- function(tab,...){ oldpar <- graphics::par(n…

Rcpp で class の array を使うメモ

コンパイル時に大量の警告が出るがとりあえず動く // [[Rcpp::plugins(cpp11)]] #include <RcppArmadillo.h> #include <array> // [[Rcpp::depends(RcppArmadillo)]] ////////// //family// ////////// class family { public: double a_gam; double b_gam; double tau_; virtual dou</array></rcpparmadillo.h>…

藤原香織、渡辺澄夫(2006)ベイズ検定による変化点検出のシミュレーション

これです: http://watanabe-www.math.dis.titech.ac.jp/users/fujiwara/doc/fujiwara_ibis2006.ppt.pdf この研究成果は論文化されているようですが, オープンアクセスではないみたいです: https://search.ieice.org/bin/summary.php?id=j91-d_4_889&catego…

変化点のあるポアソン過程へのパラメトリックな曲線あてはめ

Lambda <- function(x,a,b,c){a*x-b*log1p(exp(-x+c))+b*log1p(exp(c))} lambda <- function(x,a,b,c){a+b/(1+exp(x-c))} ll <- function(par,y,Tau){ a <- exp(par[1]) b <- exp(par[2]) c <- exp(par[3]) Lambda(Tau,a,b,c)-sum(log(lambda(y,a,b,c))) } d…

「研究仮説が正しい確率」(豊田、2017)のシミュレーション

「研究仮説が正しい確率」(豊田, 2017)の分布 - 廿TT ではやや甘の結論だったので、このエントリでは「研究仮説が正しい確率」という言葉は使うべきできない、とはっきり言い切ることを目指します。 「p 値を使って学術論文を書くのは止めよう」(https://…

多項ロジスティック回帰(ソフトマックス回帰)の変分推論をRで(Adagrad)

abrahamcow.hatenablog.com のおまけです。 AdaGradのすすめ - Qiita を参考にAdagradというのを使ってみたら、若干収束が速くなりました。 row_softmax <- function(x){ ex <- exp(x) ex/rowSums(ex) } VIlogistic <- function(Y, X, lambda=1, lr=1e-5, ma…

多項ロジスティック回帰(ソフトマックス回帰)の変分推論をRで

『ベイズ推論による機械学習入門』の例題です。 機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)作者:須山 敦志出版社/メーカー: 講談社発売日: 2017/10/21メディア: 単行本(ソフトカバー) サポートページ(GitHub - sam…

「研究仮説が正しい確率」(豊田, 2017)の分布

「研究仮説が正しい確率」が話題になっている。 https://twitter.com/genkuroki/status/1220130764148236289 一応僕の立場を明確にしておくと、事前分布を一様分布にしたとしても、事後分布の振る舞いは尤度関数の設計に依存するので、データを生成した真の…

EMアルゴリズムによる右打切りポアソン分布のパラメータ推定(Rcpp)

右打切りのポアソン分布なんてあまり例がないかもしれないが、むりやり考える。仮定1:ある機械の部品はレート のポアソンプロセスに従って繰り返し故障する。 仮定2:部品が故障したとき、ユーザーはポアソンプロセスとは独立に確率 でその機械の使用をやめ…

ベイズ的逐次ABテストの注意点

要約 ベイズ的にやったとしても、逐次的にABテストの評価をやったらアルファエラーは保たれないよ ベイズ更新 まず、ベイズ更新について簡単に説明しておきます。ベイズの定理や条件付き確率については他のページを参照してください。以降、3 つの分布が出…

Juliaの自動微分を使ってより自由なFactorization Machines

Factorization Machines の解説はこの記事がわかりやすかった: 一歩Matrix Factorization、二歩Factorization Machines、三歩Field-aware Factorization Machines…『分解、三段突き!!』 - F@N Ad-Tech Blog ただ Factorization Machines を動かすだけなら…

可視化で理解するKPIツリー:コンバージョンの寄与度分解

背景 昔Web系コンサルをやっていたときKPIツリーとかを習ったんだけどぼくはこれが苦手だった。 例えばあるサイトのコンバージョン(CV)はサイトへの訪問数(SS)×コンバージョンレート(CVR)だからコンバージョンというゴールを増やしたかったら、SSを増…

有名なポンプデータの例題(Rcppによるギブスサンプラー練習)

Pumps です。 正直これをわざわざ階層ベイズでやる意味がよくわからない。 #include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] List GibbsPomp(int iter,NumericVector x, NumericVector t, double alpha, double gamma, double delta) { int N = x.lengt</rcpp.h>…

Rcppでプロビット回帰のギブスサンプラー

なんでこんなの書いていたのか忘れてしまった。 #include <Rcpp.h> using namespace Rcpp; double rnorm_u (double eta){ return R::qnorm(R::runif(R::pnorm(0,eta,1,true,false),1),eta,1,true,false); } double rnorm_l (double eta){ return R::qnorm(R::runif(0</rcpp.h>…

r×2分割表のオッズ比をロジスティック回帰っぽく(閉形式で)計算する関数を書いた

R

R です。 使いみちがあるのかわからない。 logoddsratio <- function(mat,alpha=0.95){ cj <- log(mat[,1])-log(mat[,2]) be <- c(cj[1], cj[-1] - cj[1]) Nj <- rowSums(mat) pj <- mat[, 1]/Nj se <- pj * (1 - pj) * Nj se <- c(1/se[1], 1/se[-1] + 1/se…

比率の差を用いてA/B テストに必要なサンプルサイズをフェルミ推定

多くのABテストでは次のような分割表を評価します。 成功 失敗 a c b d ABテストを仮説検定の枠組みの中で、例数(サンプルサイズ)設計をきっちり行ってABテストを実施しようと主張する人もいます。 が、たぶん無理だと思います。 その理由として次のように…

[R + Google アナリティクス]ページ/セッションの信頼区間

ページビューをセッション数で割った値(ページ/セッション)を見ることがあります。 この数字は訪問あたりの閲覧ページ数と解釈できます。 これの信頼区間を出す方法を考えました。 ページビューをサイトからの離脱に至るまでの試行回数と考えると幾何分布…

状態空間モデルで「見せかけの回帰」を外せるか

(本文中のコードは全てRです) ランダムウォークする系列どうしで回帰を行うと全く別々に生成されたデータであっても有意差が出やすくなる。 set.seed(1) x <- cumsum(rnorm(100)) y <- cumsum(rnorm(100)) fit1 <-lm(y~x) summary(fit1) > summary(fit1) C…

ガンマ分布の形状母数が1より小さいかどうかの検定を考えた

予想 () を独立に同じガンマ分布に従う確率変数とする. とする. とする. ここで, はオイラーのガンマ定数(オイラーの定数 - Wikipedia) はトリガンマ関数 (ディガンマ関数 - 機械学習の「朱鷺の杜Wiki」) である. ガンマ分布の形状母数が1で n が十分大き…

[RStan]ディリクレ・多項モデルで μ's とAqours の人気の差を調べる

背景の整理 前に書いたエントリが引用されて嬉しかったのでもう一ネタ書きます. アイドルの人気投票における多項ロジスティック回帰の適用例 - gaiaskyの技術メモ [RStan]多項ロジスティックモデルで μ's とAqours の人気の差を調べる - 廿TT モデル キャ…

Dirichlet Multinomial モデルのトピック数を ELBO で決める

R

モデルと推定方法については以下に書いた: 変分ベイズによるトピックモデル(Dirichlet-multinomial Model)のパラメータ推定の高速化 - 廿TT ELBO の導出は後で説明します。 R によるシミュレーション library(gtools) VBDirMult <-function(Y,L=2,alpha=r…

NMF(非負値行列因子分解)の次元数を ELBO で決める

R

NMF の変分推論については非負値行列因子分解をRで(ベイズ推論による機械学習入門)2 - 廿TTに書いた。ここでは同じモデルの ELBO (evidence lower bound) を導出する。 ELBO 後で説明します。 R による実装例 NMFVB <-function(Y,L=2,alpha=1,beta=1,maxit…

GaP モデルのトピック数を ELBO で決める

R

ここでは [math/0604410] Discrete Component Analysis の Gamma-Poisson モデル(GaP)の ELBO (evidence lower bound) を導出する。行列の分解がトピックモデルの一種として解釈できることは以下に書いた: 行列の知識ゼロからはじめてトピックモデル(GaP…

ポアソン過程・フラットモデル(BTYDモデル番外編)

R

モデルと尤度 RFM 指標から将来の購買回数を予測するモデルとして、次の2つを紹介しました。 Pareto / NBD モデルの変分推論(RFM指標から将来の購買回数を予測する) - 廿TT BG / NBD モデルの変分推論(RFM指標から将来の購買回数を予測する) - 廿TT ポイ…