R

私の今後は?潜伏期・再発

曝露からの時間の関数で発病率が決まるとする 最終的に発病しないこともあるとする 曝露からある時刻まで発病していないときに、その後、発病するであろう確率を知りたいことがある いわゆる感染症で曝露がわかっているとき 癌の術後再発 遺伝性疾患で浸透率…

二項分布・ベータ分布、ポアソン分布・ガンマ分布

事後分布推定作者: ryamada発売日: 2014/08/03メディア: Kindle版この商品を含むブログ (1件) を見る 事後分布推定、プライアは、とか、色々あるけれど、ひとまず、事前情報なしでデータがあったときに、Rを叩いてみることだけを目指しての資料 Rmdファイル…

連立常微分方程式

なる常微分方程式は と分解したうえで と表される ただしは固有値について、を対角成分とする行列である こんなに単純で良ければ、が時系列データとして観察されれば、時間差分を時間微分に代用して考えて、線形回帰すればよい やってみる d <- 4 A <- matri…

医学を作るシミュレーション

医学を単純化するとどうなるかやってみる(こちらの『知識ゼロからの診療』のための模擬設定) 生物が必要 生物に多様性を持たせることが必要 病気が必要 病気の定義が必要 構造障害・機能障害(解剖・生理) 因果推論が必要 治療が必要 まずは生物の構造を設定 …

メモ

var.equal=TRUEの影響を見る ma <- 120 mb <- 120 sa <- 10 sb <- 10 + (0:(-9)) na <- 15 nb <- 10:15 sb.nb <- expand.grid(sb,nb) n.iter <- 10000 ps <- matrix(0,length(sb.nb[,1]),n.iter) for(i in 1:length(sb.nb[,1])){ for(j in 1:n.iter){ x <- r…

ベータ分布の交叉点

ベータ分布が二つあるとする その累積分布を考える。ベータ分布の累積分布は、regularized incomplete beta 関数 2つのベータ分布の累積分布関数の交点を計算してみる # x1は一つ目のベータ分布を決める、二つの成否回数ベクトル # x2は二つ目のベータ分布…

即効データシミュレーション実習

Rの初心者だけれど、「もう大人」なので、なんとか使えるようになろうという意識は高い人には、「つまらない」例ではなくて、少々わかりにくくても、「できるといいな」という例がよい こちらのJuly15-の週の" Click to download today's text File " こんな…

Dirichlet diffusion trees

先日の記事の続き 昨日の木はDirichlet diffusion treeの分岐プロセスをグラフオブジェクトに乗せて作る話だった 昨日の記事では、分岐確率を時刻の関数とし、さらに、各エッジで同道している「先行者」の数で重みづけすることを述べたが実装は面倒なので省…

Dirichlet diffusion trees

先日の記事の続き 作り方 原点から第一の要素がランダムウォークする。単位時間後とすれば、ランダムウォークするので、平均0、分散がある値(1としてもよいだろう)の点(t=1,x)に到達する 第二以降の要素は、原点から出発し、「分離」イベントが起きるまでは…

Dirichlet diffusion trees

資料1 資料2 資料3 資料4 データがあって、このモデルに合うか、あのモデルに合うか、この仮説に合うか、あの仮説に合うか、というのを考えるとき「尤度」「尤度比」を使って、事前確率を事後確率に変える このモデル・あのモデルと考えつつ、MCMCで事前…

Latent class mixed effect model

参考はこちら library('lcmm'); data(data_hlme) head(data_hlme) データは経時的に観察するYと観察時刻Timeが各行で異なり、それ以外の共変量は個人の番号IDと同様観測回数ごとに繰り返して入力されている形式。以下では、ID=(1,2)の2人分のデータが見えて…

0の行・列を縮める

複数のカテゴリカル変数が作る、多元表があって、各変数の周辺度数が0になっているかもしれない(のだが、元の多元表の形で考えていたい)とき、表を縮めたり(そのあと、縮めることで処理が容易になることをしたり)、いったん縮めた表を元の形に戻したりしたい…

Quantile Regression

Wikiはこちら 日本語記事で扱っているのはこちら 観察データがあったとする クオンタイルの回帰(直)線を引く前に、クオンタイル値の推定について考える 今、-クオンタイル値をとしてうまいこと求めたい が-クオンタイル値であるとき、観測値がより小さければ…

Topology Data Analysis

PDF 以前もメモしていた(こちら)

少し改良

R

カテゴリ数を増やしつつ、その「未観測カテゴリ」の確率をプールして計算するだけなら をカテゴリ数仮説の相対尤度とし、「未観測カテゴリ1個分」の生起確率期待値はなので、の未観測カテゴリ数分を合算すればについて、カテゴリ数仮説の相対尤度重みづけ平…

次に何が出る? ディリクレ、色々お試し編

一昨日・昨日と、カテゴリ数が未知のときのディリクレ分布の合算の話をしてきた( まだ、この方法でOK、という確信はないのだが、まあ、これまでのところ、出てくる値もそれなりに妥当そうなので、このまま行ってみる(一番、危なそうなのは、事前確率をすべ…

次に何が出る? ディリクレ実践編

本当にできるかやってみよう まずは、kカテゴリ観察されているときに、元集団のカテゴリ数がkである、と限定した場合 この場合は、kカテゴリでのディリクレ分布そのものであり、ディリクレ分布に従うとき、n回のうち回観察されたカテゴリの期待値はであると…

次に何が出る? ディリクレお話編

昨日は、複数のカテゴリからなる「元集団」から有限回サンプリングしたら、k種類のカテゴリが観測された場合の話をした。 そのときに観測されたk種類のカテゴリだけが元集団を構成しているとみなして、元集団のカテゴリ割合をディリクレ分布として推定する場…

ディリクレ分布をいじる

kカテゴリあって、全部でN回サンプリングしたところ、それぞれが回ずつ観測されたとする。ただし このとき、サンプリング元のkカテゴリの割合をは確率変数であって、そのような確率はに従うと想定するとよいことが知られていて、この分布をディリクレ分布と…

Ewen’s sampling formulaとクラスタ数

http://www.stat.uchicago.edu/~pmcc/pubs/clusters.pdf クラスタ数が不明なときに、観測データからクラスタ数をどう読み取るか Ewens's sampling formula(Wiki) 不定のクラスタをもつ集団からの有限サンプリングに関する集団遺伝学・生態学的な式 untbパッ…

離散ラプラス分布

Y染色体は染色体全体のハプロタイプを問題にする、という特別な事情があるために、集団にどんなハプロタイプがどれくらいの割合であるのか、とか、二つのY染色体ハプロタイプがたまたま一致する確率はどれくらいか、と言うことが問題になります。 しかしなが…

リプレース

メモリンク n <- 50 n.iter <- 100000 N <- 40 r.with.replace <- r.without.replace <- matrix(0,n.iter,N) for(i in 1:n.iter){ x <- rnorm(n) r.with.replace[i,] <- sample(x,N,replace=TRUE) r.without.replace[i,] <- sample(x,N,replace=FALSE) } m.w…

Non-negative Matrix Factorization

こちらの論文で腫瘍のシークエンスデータを用いて、変異の入った遺伝子とそうでない遺伝子との情報を行列化し、それにNon-negative Matrix Factorizationを使って腫瘍のクラスタリング、特性(治療反応性とか)予測、特性に影響しているパスウェイ検出などを行…

順序とかのための正確検定法パッケージ

量的データだけれど、分布を仮定していいのかわからない 同じ値がある(Tie(タイ)がある) グレード別のデータ(タイばっかり) というような場合 しかもサンプル数が、え、これだけ? いくつかの方法を正確検定中心にまとめたらしいパッケージがexactRankTests…

Association Ruleの視覚化パッケージarulesViz

arulesパッケージとarulesVizパッケージとを使ってAssociation Rule法の出力の中身とその視覚化について確認します。 以下はRmdファイルです。こちらなどを参考にすれば、自由にhtml文書化・epub化できます(が、それが面倒な場合には、kindleで→こちら) メモ…

rstan

Stanについてはこちらに紹介記事があり、環境整備に関しては、こちらを参考にして、以下のようにやればよいようです(Windowsの場合) 基本的には Rがあって Rtoolsを入れて、Rcppパッケージが動くようにして、inlineでのRcppによるRのC++化を可能にし rstanを…

対応の壊れたゲーム結果を復元する

R

12人の参加したボウリング大会 全員が4ゲームを実施して記録した 第1ゲームでの投球順の記録はある 第2ゲーム以降では第i投球者がどのようなスコアを残したかはわかるが、各ゲームで誰が何番目に投げたかがわからなくなった とはいえ、プレーヤごとにかな…

二次形式なLoss fucntion

こちらでGaussian Sequence Modelをなぞっている 遅々として進まないのだが、まあそれはよいとして、今は、Loss functionとして二次形式なもので正定値な行列を挟んだものを使うと、Bayesian estimatorは観測値に基づく真値の事後分布の期待値になるよ、それ…

前期 医科学のためのデータシミュレーション

R

こちらの続き 医学・生物学データの扱いに慣れるための、データシミュレーションセミナーの資料(医学研究科大学院コース ゲノム遺伝学@京大 2014:こちら) htmlファイルはこちら(Login as a guest でログインすると見られます) こちらを参考にすればhtml化可…

医科学のためのデータシミュレーション 3

R

こちらの続き 医学・生物学データの扱いに慣れるための、データシミュレーションセミナーの資料(医学研究科大学院コース ゲノム遺伝学@京大 2014:こちら) htmlファイルはこちら(Login as a guest でログインすると見られます) こちらを参考にすればhtml化可…