R

ディリクレ過程混交正規分布様球面分布

球面に粗密がランダムな点を取りたい ディリクレ過程で定まる無限項の多項分布を想定し、個々の項に対応して球面正規分布様のそれを取る ただし、球面積分布様のそれ、とは、球面上の一様乱点を接点とする接面に、接点を中心とした2次元正規分布をとる。その…

疑似1細胞発現データを作る

複数の細胞からのデータがある 細胞数が分かっている リード数を細胞数に分けて、疑似1細胞発現データを作りたい sample()関数を使って疑似作成できることを示し # 検体数 N # ある遺伝子のNGSリード数が Mだった # m1 + ... + mN = M; mi >= 0となるような…

Traceが1の自己随伴行列。そのpure stateとmixed state

Traceが1の自己随伴行列で固有値がすべて非負のものは、密度行列と呼ばれ、固有値が確率を、固有ベクトルが対応する(離散的)状態を表している このような行列のうち、単位ベクトルの複素共役クロネッカー積([tex:|x> Methods of Information Geometry (Ta…

正準相関解析、stats::cancor()とcandisc::cancor()

解説記事

医療経済分析

Rに個人の時系列変化を追いかけて治療法(治療薬)の経済性を解析するパッケージ hesim がある。そのVignette(こちら)のコードをなぞってみる。オリジナルサイトのコードだとエラーが出るので、ちょっと修正してある

パワー計算の実習

ランダム家系データ作成

対立仮説の下でのSNP table

t-SNEのt-分布

次元削減の1手法のtSNE 高次元の標本間の遠近情報を ある点を中心として別の点が観察される確率として評価する その確率を、観測標本全体で重み付き標準化する さらに、その標準化確率行列が非対称なので対称化する 低次元に埋め込むとき 全く同様に低次元…

デンドログラムを二次元平面の上に描く

二次元平面上に標本があり、それらの間に階層的クラスタリング構造が得られたとする 二次元平面上に標本を配置しつつ、その平面から立ち上がるようにデンドログラムを描きたいかもしれない library(FactoMineR) # Compute PCA with ncp = 3 res.pca <- PCA(U…

球面調和関数係数も一般化逆行列で推定

球面に場があったとする それを球面調和関数で分解する 球面調和関数は直交基底関数なので、個々の関数が場を説明する成分を積分して計算してもよいが、SPHARMというツール(こちら)では、一般化逆行列を使っている(この論文のメソッドセクション) どういうこ…

ノンパラメトリック・ベイズ実践編2 回帰

まず、単純に線形回帰 s <- rnorm(n,0,1) head(s) m <- matrix(s,ncol=2) head(m) lm_res <- lm(m[,2] ~ m[,1]) lm_res plot(m) abline(lm_res) pythonで行こう Scikit learn の gaussian process regression そのExamplesの一つ カーネル関数から、値のsimi…

ノンパラメトリック・ベイズ実践編

Non-parametric bayesian clustering Data set simulation n <- 1000 d <- 4 k <- 5 s <- sample(1:k,n,replace=TRUE) m <- matrix(rnorm(d*k,sd=6),k,d) X <- matrix(0,n,d) for(i in 1:k){ tmp <- which(s == i) r <- matrix(rnorm(d * length(tmp))) X[tm…

固有値分解する

の最小化 の最小化 Mの最小固有値に対応する固有ベクトルが一番拡大率が小さいから、その固有ベクトル方向のが求める解 n <- 100 X <- matrix(rnorm(n*2),ncol=2) Y <- X %*% c(1,2) + rnorm(n,0,0.01) library(rgl) plot3d(X[,1],X[,2],Y) lm(Y ~ X-1) XY <…

Maximum Mean DIscrepancy その3

Rのkernlab library(kernlab) # create data x <- matrix(runif(300),100) y <- matrix(runif(300)+1,100) mmdo <- kmmd(x, y) mmdo > mmdo Kernel Maximum Mean Discrepancy object of class "kmmd" Gaussian Radial Basis kernel function. Hyperparameter…

特性関数

確率変数Xがあるという。たとえば正規分布に従う変数。 確率密度関数が書けたりする 今、Xと関係する別の確率変数Yを考える ただし、Yは 複素数である 実数変数tによって変わるものとする 実際と定める この複素確率変数には「平均〜期待値」がある この期待…

kernlab パッケージ

カーネル法の代表格がサポートベクターマシン サポートベクターマシンの実装はC/C++,MATLABにもあり、Rでもe1071にそのカウンターパートがある サポートベクターマシン以外にもカーネル法の使い道はある カーネル法一般利用をするための諸関数を提供するのが…

MMDでChIP seq解析

MMDiff2パッケージ source("https://bioconductor.org/biocLite.R") biocLite("MMDiff2") library('MMDiff2') library('MMDiffBamSubset') ExperimentData <- list(genome='BSgenome.Mmusculus.UCSC.mm9', dataDir=system.file("extdata", package="MMDiffBa…

log sum exp

メモ 指数関数を使うと値が小さくなりすぎるものの足し合わせの際に、複数の値の相対値がそれほど小さくならないことを利用して、計算の精度を担保する方法 f <- function(a,b){ #log(exp(a)+exp(b)) b + log(exp(a-b)+1) } g <- function(as){ if(length(as…

Mendelian Randomization

以前、こんな風にまとめた 別の書き方をしてみる 簡単に書き直すと、YにXが影響を与えているのだが、それは隠れた因子UがXとYに影響を与えているからかもしれないので、Uに関係なくXがYに影響を与えているかを知りたい。しかしながらUはどんなものがあるのか…

ぱらぱらめくる『DPpackage: Semi- and Non-parametric Modeling in R』

文書 1. Introduction Semiparametric Bayes と Nonparametric Bayes (BSPとBNP) Bayesian アプリ。RでDirichlet過程を実装しているのは(DPpackageの他に)bayesmがあるが、扱えるモデルは限定的 DPpackageが実装しているモデル DP mixtures of DP (Antoniak'…

ディリクレ過程で二項分布 その2

昨日の記事でRのDPpackageのDPbetabinom()関数の使い方をさらった 観測データが320のオブジェクトについて9回繰り返し試行をしているわけだが、9回の繰り返し試行で、それの「表確率」の推定はそれほど精度がよくできるものではない。 個々のオブジェクトで…

ディリクレ過程で二項分布

一昨日の記事でRのDPpackageに触れた 今日の記事では、二項分布・ベータ分布を例にいじってみる 関数はDPbetabinom()である まずは、この関数のヘルプ記事のExamplesから 全部コピペして動かしてみる library(DPpackage) # Data data(rolling) y <- cbind(ro…

DPpackage

メモ メモ2 DPDensity() 関数の場合 正規分布の混合を考える 構成正規分布数を不定にするところがノンパラメトリック i番目の正規分布: : 各正規分布を定めるパラメタは、ディリクレプロセスで発生させる。ただし、そのディリクレプロセスはで決まる。G0か…

Gaussian Process Regression

こちらにGaussian Process Regressionをベタコードで解説してある 手元のR環境だと、描図エラーが出るので、エラーを出す行だけコメントアウトしたものを以下に再掲する xの座標ペアごとにどれくらいの関連を入れるかを分散共分散計算関数で指定し 観測デー…

色覚異常フレンドリーなカラーパレット

色覚異常には先天性のものと後天性のものがあるが、先天性のものにも複数の種類が知られている。日本人に多いのは先天性赤緑色覚異常で、2種類に分かれます(こちら)。英語での名称はDichromacy (2種類がprotanopia and deuteranopia) 色覚異常がどういうこ…

rstan

ようやくrstanがわかってきたのでまとめ直す まとめ データがある 確率モデルを定め、対数尤度計算をコード化する rのc++連携の仕組みを使ったrstanパッケージのstan()関数を使ってStan推定をする 何がよいのか 上記のまとめは読んでもつまらないまとめにな…

RとSTANとでBayesianの基礎

ベイズ推定をする データがある モデルがある モデルは、パラメタを持っていて、そのパラメタの値を定めると、データを観察する確率が定義できる Stanでは、このパラメタがどういう分布を取っているかを乱数を使って標本分布として返す パラメタ以外にも、ど…

rstanの調べもの

modelの部分が何をやっているのか…対数尤度を計算している。その書き方の基礎の基礎→こちら stanの基本的なことだけれど、なかなか探すのが難しいことをかいつまんで書いてあるサイト→こちら 一般化線形回帰を使ってrstanの動かし方を確認する→こちら stanフ…

Rで探し物

R

21の探し物ツール Navigation; R-bloggers; Task views; Rdocumentation.org; sos package; Rコマンド・関数 ??; apropos; ls; methods; getAnywhere; 隠れた関数も表示 :::; find; args; grep; %in%; str; getwd; file.choose; Spyglass summary; browser; …