R

Metropolis Hastings サンプリングする

# 値の台 xx <- seq(from=-5,to=20,length=1000) # 適当な(密度関数) dd <- dnorm(xx) + dnorm(xx,1,2) + dnorm(xx,10,2) # 確率密度関数っぽい形 plot(xx,dd) # library(mcmc) # 確率密度に比例したあたいの対数を返す関数を書く h <- function(x){ log(dno…

グレブナー基底で連立多項式を解く、量子計算機でグレブナー基底を探す

https://arxiv.org/pdf/1903.08270.pdf pythonのsympyを使うとグレブナー基底が計算できる sagemath on pythonでもできる from sympy import * x,y,z = symbols('x y z') eqs =[x**2 + y**2 + z**2 - 4, x**2 + 2 * y**2 - 5, x*z-1] result = groebner(eqs,…

numpy.ndarray と Rのapply()

Rの行列で、各行の和を出そうと思ったら M <- matrix(rnorm(100),10,10) apply(M,1,sum) だろう 同じことをPythonでやるのはどうするのかな、と思った Pythonで行列を扱うならnumpyがよいのでそうすることにしたとして 値のセットの和を出すsum()関数を、軸…

RとPythonの両方で使えるplotlyでの3Dスキャタープロットを覚える

Rもpythonもjupyter notebookで文書作成できる 3次元プロットをしなくてよいのなら、Rでは普通の二次元プロッティング関数を使えばよいし、Pythonではmatplotlibを使えばよい 3次元プロットを描いてぐりぐり動かすとなると、Rではplot3dを使い、pythonでは…

混合正規分布のEMアルゴリズム推定

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

球面に粗密がランダムな点を取りたい ディリクレ過程で定まる無限項の多項分布を想定し、個々の項に対応して球面正規分布様のそれを取る ただし、球面積分布様のそれ、とは、球面上の一様乱点を接点とする接面に、接点を中心とした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の座標ペアごとにどれくらいの関連を入れるかを分散共分散計算関数で指定し 観測デー…