ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2016-06-12 ぱらぱらめくる『統計データ分析 ―ベイズ的〈ポストp値時代〉の統計

[][][]準備:ぱらぱらめくる『統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―』

  • 準備
  • こちらに書籍用の諸ファイルがあります
  • ダウンロードしたzipファイルを展開すると、
    • 3つのフォルダ
      • myfunc
      • srcA
      • stan
    • 2つのファイル
      • Readme.txt
      • RstanStan.pdf
  • があります
  • Readme.txtは、このフォルダとその中の諸ファイルを使うための準備について書いてあります
    • 必要なものは
      • R本体
      • rstanというパッケージ
install.packages("rstan")
        • でインストールします
      • Rtools(という道具一式。c,c++のコンパイルなどをしてくれます)
  • その準備が揃ったところで
  • Working directoryを、このzipを展開したところ(statという名前のディレクトリ)に設定します
    • そのstatディレクトリを指定します。Desktopに展開したなら、たとえば
setwd("C:/Users/ryamada/Desktop/stat")
    • のように。
  • Rstanのインストールについてはこんなページもあります

[][][]実行:ぱらぱらめくる『統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―』

  • 実行は、srcAディレクトリのファイルをコピーペーストして実行するだけ
  • 6章分
source('myfunc/myfunc.R')
  • 冒頭でmyfuncディレクトリのmyfunc.Rソースを読み込みます
  • このmyfunc.Rは、myfuncディレクトリ内のそのほかのいろいろなxxx.Rソースを読み込みます
  • 第1章用のファイル"part1.R"はstanを使うより前
  • 第2章用から、stanを使います。それは第2章用ファイル"part2.R"の冒頭に
library(rstan)
  • としていることからわかります
  • どこでstanをしているかと言えば、part2.Rの中で使われている
#正規分布に関する推測
out <-G1mean(x,prior=T,mL=0, mH=1000, sL=0, sH=100)       
  • にあるG1mean()という関数が見慣れません
  • これはmyfuncディレクトリ内にある関数で、以下のように中でrstanパッケージを呼び出したりしており、stanによる推定処理を決めています
##########################################################################
#1群の正規分布に関する推測
#■入力
#x:ベクトル形式のデータ 
#prior=F : 論理値。Tなら事前分布の範囲を指定。Fなら指定せず。
#mL, mH, sL, sH : 事前分布のパラメタ
#prob=c(0.025, 0.05, 0.5, 0.95, 0.975): 事後分布で報告する確率点
#see=1234,cha=5,war=1000,ite=21000 : MCMC関連のパラメタ
#■出力
#fit:stanの出力, par:母数リスト, prob:確率ベクトル, 
#mu:平均, sigma:標準偏差, xaste:予測分布, log_lik:対数尤度
##########################################################################
G1mean <- function(x,prior=F, mL=-1000, mH=1000, sL=0, sH=100,
prob=c(0.025, 0.05, 0.5, 0.95, 0.975),see=1234,cha=5,war=1000,ite=21000)
{
   library(rstan)
   if (prior) {
     dat <- list(n=length(x),x=x,mL=mL,mH=mH,sL=sL,sH=sH)
     scr <- "stan/G1meanPT.stan"                           # Stan ファイル名
   } else {
     dat <- list(n=length(x),x=x)
     scr <- "stan/G1meanPF.stan"                           # Stan ファイル名
   }
   par<-c("mu","sigma","xaste","log_lik")        # sampling変数
   fit<-stan(file=scr,data=dat,pars=par,seed=see,chains=cha,warmup=war,iter=ite)
   ext<-extract(fit, c("mu","sigma","xaste","log_lik"))
   mu<-ext$mu; sigma<-ext$sigma; xaste<-ext$xaste; log_lik<-ext$log_lik
   out<-list(fit=fit,par=par,prob=prob,mu=mu,sigma=sigma,xaste=xaste,log_lik=log_lik)
   class(out)<-'G1mean'
   return(invisible(out))
}
  • その内部で
scr <- "stan/G1meanPT.stan"                           # Stan ファイル名
  • という行がありますが、stanディレクトリ内の"G1meanPT.stan"というファイル(Stanのやり方を指定するファイル)を読み込んでいます
  • この.stanファイルは、以下のようになっていて、ある書式に従って書かれているらしきことが読み取れます。
data {
int<lower=0> n;                                  #データ数
real<lower=0> x[n];                              #データ
real mL; real mH; real sL; real sH;  #事前分布
}
parameters {
real<lower=mL,upper=mH> mu;             #平均(範囲指定)
real<lower=sL,upper=sH> sigma;          #標準偏差(範囲指定)
}
transformed parameters {
}
model {
x ~ normal(mu,sigma);                   #正規分布
}
generated quantities{
real xaste;
real log_lik;
xaste<- normal_rng(mu,sigma);           #予測分布
log_lik<-normal_log(x,mu,sigma);        #対数確率
}

  • 実際、この.stanファイルを引数として受け取っているのは、G1mean()関数の中の次の行
fit<-stan(file=scr,data=dat,pars=par,seed=see,chains=cha,warmup=war,iter=ite)
  • で、これによって、推定結果がfitに返っています

2016-05-21 2のd乗組み合わせデータ構造のための情報幾何

[][][][][][]Dually Flat Manifolds

  • 昨日、階層構造を持つ確率分布のための情報幾何について少し書いた
  • 特に、¥{0,1¥}^dなd元分割表のためのe-平坦、m-平坦なパラメタの取り方について、丁寧に確認してみたい
  • Rでパラメタ変換をしてみるのが、手っ取り早そうなので、それでやってみる
  • Rmd
---
title: "Dual flat parameterization for {0,1}^d"
author: "ryamada"
date: "Friday, May 20, 2016"
output: html_document
---

# Brief introduction

{0,1}の値を取るd個の変数$X_i;i=1,2,...$があるときに、その値セットの取り方($\mathbf{x} = (0,0,...,0),(1,0,...),...,(1,1,...,1)$)は
$D=2^d$通りある。

今、この$D$通りの取り方の生起確率ベクトルを$freq=(f_1,...,f_D)$とすれば、
$\sum f_i = 1$の制約があり、自由度は$D-1$である。


このような階層構造を持つ組み合わせについては、情報幾何的な観点から有用な、パラメタの取り方が2通り知られている。

It is know that there are tow useful ways to parameterize the frequency with the hierarchical structure.

この文書の目的は、頻度ベクトルを2通りのパラメタに置き換えたり、パラメタの値を指定して、頻度ベクトルを作ったりする道具を提供すること。

This document is written to offer utility functions to transform frequency vectors to two types of parameters and their inverse.




詳しい理論は ["Information Geometry on Hierarchy of Probability
Distributions" (IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 5, JULY 2001)](https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwja3aaXwOXMAhXHXaYKHVSXCzAQFgghMAA&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel5%2F18%2F20133%2F00930911.pdf%3Farnumber%3D930911&usg=AFQjCNFd9avfcHWB3jmeVpfvnuWzCKo5rQ&sig2=OBOG6MXZgpPJvHd2JksQIQ) を参照のこと。

You can learn its theory in ["Information Geometry on Hierarchy of Probability
Distributions" (IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 5, JULY 2001)](https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwja3aaXwOXMAhXHXaYKHVSXCzAQFgghMAA&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel5%2F18%2F20133%2F00930911.pdf%3Farnumber%3D930911&usg=AFQjCNFd9avfcHWB3jmeVpfvnuWzCKo5rQ&sig2=OBOG6MXZgpPJvHd2JksQIQ)

## e-flat parameterization

$\mathbf{\theta}$は長さ$D-1$のパラメタベクトル。

$\mathbf{\theta}$ is a parameter vector with $D-1$ elements.

$\mathbf{\theta} = ((\theta_1,\theta_2,...,\theta_d,),(\theta_{1,2},...,\theta_{d-1,d}),...,(\theta_{1,2,...,d}))$

$\mathbf{x}=(x_1,...,x_d)$というパターンを取る確率
$$
\log{p(\mathbf{x})}=(x_1,...,x_d))= \sum \theta_i x_i + \sum \theta_{i,j} x_i x_j + ... + \theta_{1,2,...,d} x_1 x_2 ... x_d - \psi
$$

$psi$の値は$\mathbf{\theta}$が決まると(自由度がもうないので)一意に決まる。

The value of $\psi$ is not free because $\mathbf{\theta}$ has $D-1=df$ elements.

$\theta$ は1要素の寄与、要素ペアの寄与、要素トリオの寄与…に対応づく。


$\theta$ represents single effect, pairwise effect, triowise... and so on.

## m-flat parameterization

$\mathbf{\eta}$は、長さ$D-1$のパラメタベクトル。

$\mathbf{\eta}$ is a parameter vector with $D-1$ elements.

$\mathbf{\eta} = ((\eta_1,\eta_2,...,\eta_d,),(\eta_{1,2},...,\eta_{d-1,d}),...,(\eta_{1,2,...,d}))$

$\eta$の各要素は、変数(単独・組み合わせの)期待値を表す。(多元分割表で言えば)周辺頻度にも対応する。

$\eta$ corresponds to the expected value of single/combination of variables, or in other words, marginal frequency of multi-way table.

$$
\eta_{i} = Prob(x_i=1)\\
\eta_{i,j} = Prob(x_i=x_j=1)\\
...
\eta_{1,2,...,d} = Prob(x_1=x_2=...=x_d=1)
$$

# R codes in the ryamada22/Ronlyryamada github-package

# Utility functions and their usage example
```{r}
#library(devtools)
#install_github("ryamada22/Ronlyryamada")
library(Ronlyryamada)
example(my.haplotypes)
```

# Vignettes on functions

## Enumeration of $2^d$ patterns
d個の変数、$2^d$の場合の列挙。
1行目から順番に、値が1である変数の数が、0,1,...,dと増えるような並び方になっている。

$2^d$ patterns for d variables.
From the top row to the bottom, number of variables with value 1 increases like 0,1,...,d.


```{}
my.haplotypes <- function(d){
  n <- 2^d
	X <- matrix(0,d,n)
	cnt <- 1
	for(i in 1:d){
		tmp <- combn(d,i)
		for(j in 1:length(tmp[1,])){
			X[tmp[,j],cnt+j] <- 1
		}
		cnt <- cnt + length(tmp[1,])
	}
	return(X)
}
```
```{r}
d <- 3
my.haplotypes(d)
```

## Matrices to convert between frequency, theta and eta.

```{}
my.freq2eta.mat <- function(d){
  X <- my.haplotypes(d)
	n <- 2^d
	M <- matrix(0,n,n)

	for(i in 1:n){
		M[,i] <- apply(X >= X[,i],2,prod)
	}
	return(t(M[2:n,2:n]))
}
my.theta2freq.mat <- function(d){
  K <- my.freq2eta.mat(d)
	n <- 2^d
	M <- matrix(0,n,n)
	M[2:n,2:n] <- K
	M[,1] <- -1
	return(M)
}
my.freq2theta.mat <- function(d){
  M <- my.theta2freq.mat(d)
	return(solve(M))
}
# eta2freq.mat = inverse of freq2eta.mat
```
```{r}
my.freq2eta.mat(d)
my.theta2freq.mat(d)
my.freq2theta.mat(d)
```

## Transformations among three parameterizations

"freq", "$\theta$", "$eta$"3通りのパラメタ表現。

"$\theta$"の場合、自動的に決まるD番目の値を計算しておく方がよい。ほかの2つのパラメタシステムの場合は足して1になるところからの差、と言う形で簡単に計算できるが、$psi$はそうではなく面倒くさいから。


"freq" stands for regular frequency vector with $2^d=D$ elements.

"$theta$" stands for e-flat parameters.

"$eta$" stands for m-flat parameters.

For $\theta$, the residual parameter "\psi" should be also calculated because it is a bit difficult than other two systems' D-th value, even though it is automatically determined. 

```{}
my.theta2psi <- function(theta){
  len <- length(theta)
	d <- log(len+1,2)
	M <- my.theta2freq.mat(d)
	#tmp <- M %*% c(0,theta)
	tmp <- my.Inf.matrix.mult(M,c(0,theta))
	delta <- sum(exp(tmp))
	psi <- log(delta)
	return(list(psi=psi,M=M))
}
#' @export

my.theta2freq <- function(theta){
	out <- my.theta2psi(theta)
	psi.theta <- c(out$psi,theta)
	M <- out$M
	return(exp(my.Inf.matrix.mult(M,psi.theta)))
}
#' @export

my.freq2thetapsi <- function(freq){
	d <- log(length(freq),2)
	M <- my.freq2theta.mat(d)
	#return(M %*% log(freq))
	return(my.Inf.matrix.mult(M,log(freq)))
}
#' @export

my.freq2eta <- function(freq){
	d <- log(length(freq),2)
	K <- my.freq2eta.mat(d)
	return(K %*% freq[-1])
}
#' @export

my.eta2freq <- function(eta){
	d <- log(length(eta)+1,2)
	K.inv <- solve(my.freq2eta.mat(d))
	tmp <- K.inv %*% eta
	return(c(1-sum(tmp),tmp))
}
#' @export

my.theta2eta <- function(theta){
	freq <- my.theta2freq(theta)
	return(my.freq2eta(freq))
}
#' @export

my.eta2thetapsi <- function(eta){
	freq <- my.eta2freq(eta)
	return(my.freq2thetapsi(freq))
}
```

## Infinite values の扱い

freqベクトルに0が含まれると、$\theta$が無限大になる。

無限大を含めて計算できるようにするために、少し、工夫が必要。そのユーティリティ関数。

無限大と非無限大とを分離して行列計算し、最後にまとめる。

```{r}
my.matrix.Inf.separator <- function(M){
  if(!is.matrix(M)){
		M <- matrix(M,ncol=1)
	}
	M.inf <- M.noninf <- matrix(0,length(M[,1]),length(M[1,]))
	infs <- which(is.infinite(M))
	noninfs <- which(!is.infinite(M))
	M.noninf[noninfs] <- M[noninfs]
	M.inf[infs] <- sign(M[infs])
	return(list(M.noninf=M.noninf,M.inf=M.inf))
}

my.Inf.matrix.mult <- function(M1,M2){
	out1 <- my.matrix.Inf.separator(M1)
	out2 <- my.matrix.Inf.separator(M2)
	
	tmp.n.n <- out1$M.noninf %*% out2$M.noninf
	tmp.i.i <- out1$M.inf %*% out2$M.inf
	tmp.n.i <- out1$M.noninf %*% out2$M.inf
	tmp.i.n <- out1$M.inf %*% out2$M.noninf
	
	tmp.any.i <- sign(tmp.i.i + tmp.n.i + tmp.i.n)
	non.zero <- which(tmp.any.i!=0)
	tmp.any.i[non.zero] <- tmp.any.i[non.zero] * Inf
	tmp <- tmp.n.n + tmp.any.i
	return(tmp)
}
```

2016-05-20 Dually Flat Manifolds

[][][][][][]Dually Flat Manifolds

  • 資料
  • 確率分布を点としてもつ多様体を統計多様体という
  • 確率分布をn個のパラメタで表すことにすると、その統計多様体はn次元多様体
  • 多様体には多様体固有の特徴があって、それはリーマン計量だったりするわけで、局所座標系を入れる・入れないに関わらずその多様体の広がり具合というのは決まっている
  • 統計多様体では、座標系を入れる。パラメタで分布を表現することと同じ
  • パラメタを入れて、多様体全体をパラメタを使った座標表示にすると、格子が張り付いたことになる
  • この格子は、うまく張り付けると、ねじれたり捩れたりしないで素直に多様体上を流れていく
  • そんなよい感じの座標表示を"flat"と呼ぶ
    • 具体的には、3次元以上の多様体のときに、任意の3パラメタの間に延び広がり方・曲り方に偏りがないようなもののこと
    • 2通りのflatがあって、e-flat, m-flatと言うものがあり、この2つのflatは相互に密接な関係がある
    • またe-,m-flatsを表す式もよく似ている
    • e-flatの定義:確率分布が¥mathbf{¥theta}で表されている
      • E¥[¥frac{¥partial^2}{¥partial ¥theta_i ¥partial ¥theta_j}¥log{p(x,¥mathbf{¥theta})} ¥frac{¥partial}{¥partial ¥theta_k}¥log{p(x,¥mathbf{¥theta})}¥]=0
    • m-flatの定義:確率分布が¥mathbf{¥eta}で表されている
      • E¥[¥frac{1}{p(x,¥mathbf{¥eta})}¥frac{¥partial^2}{¥partial ¥eta_i ¥partial ¥eta_j}p(x,¥mathbf{¥theta}) ¥frac{¥partial}{¥partial ¥eta_k}¥log{p(x,¥mathbf{¥eta})}¥]=0
  • e-flatならm-flat、m-flatならe-flat
  • e-,f-両方のflatなパラメタ表現があるとよい(のだが、そうそう簡単に見つからないかもしれない)
  • ありがたいことに…
    • 多くの確率分布のパラメタ表現は、e-flatになるような表現で書きだせる。また、そのような表現を持つ確率分布を指数分布族と呼ぶ。そのような表現は
      • p(x,¥mathbf{¥theta})=e^{¥sum ¥theta_i k_i(x)- ¥psi (¥mathbf{¥theta})}
        • ただし、e-flatな分布関数が、すべて指数分布族なわけではない
    • また、m-flatになるような確率分布族も知られている。混合分布族
      • p(x,¥mathbf{¥eta})=¥sum ¥eta_i q_i(x) + (1-¥sum ¥eta_i) q_0(x)
        • ただし、m-flatな関数が、すべて混合分布族なわけでもない
  • e-flatなパラメタ表現がわかっても、m-flat表現がわからないとつまらないし、その逆も然り
    • ありがたいことに、ルジャンドル変換と呼ばれるルールで見つけられる。以下の変換ルールを見ればわかるように、¥mathbf{¥theta}¥mathbf{¥eta}とは、要素ごとに対応関係がある
      • ¥eta_i = ¥frac{¥partial}{¥partial ¥theta_i} ¥psi(¥mathbf{¥theta})= E¥[k_i(x)¥]
      • ¥theta_i = ¥frac{¥partial}{¥partial ¥eta_i} ¥phi(¥mathbf{¥eta})
        • ¥phi(¥mathbf{¥eta}) = E¥[¥log{p(x,¥mathbf{¥eta})}¥]
  • 対数尤度関数での関係
    • 確率分布を2通りのflatなパラメタ表現をすること、それらの相互関係について書いてきた
    • 確率分布関数を尤度関数とみなすとき、尤度関数の代わりに対数尤度関数を用いることも多い。確かに、上記の記載でもあちこちで対数化した関数が用いられていた
    • 以下では、対数尤度関数¥log{p(x,¥mathbf{¥theta})},¥log{p(x,¥mathbf{¥eta})}としてやり、そのパラメタ方向偏微分を成分とするベクトルを考えると、¥theta,¥etaとが「正規直交基底関係」のように見える(実際は、¥mathbf{¥theta}のそれぞれの成分同士は互いに直交ではなく、¥mathbf{¥eta}も同様なので、e-flat側とm-flat側との対応がきれいだ、ということになる
    • <e_i,e^*_j> = ¥delta_{ij}
  • Dual flatなパラメタ表記を用いると、次が言える
    • ¥psi(¥mathbf{¥theta}) + ¥phi(¥mathbf{¥eta}) + ¥mathbf{¥theta} ¥cdot ¥mathbf{¥eta} = 0
    • この式は、p(x,¥mathbf{¥theta}),p(x,¥mathbf{¥eta}という、「同じ関数」「多様体上の同じ点」の異なるパラメタ表記に関すること
    • ここで、「異なる点」についてこれを用いる
    • D¥[p:p’¥]=¥psi(¥mathbf{¥theta}) + ¥phi(¥mathbf{¥eta}’) + ¥mathbf{¥theta} ¥cdot ¥mathbf{¥eta}’というものが定義できて、これはp=p’のときに0でそうでないときには正の値をとる
    • これがKL divergenceとなっている
    • D¥[p:p’¥] = E_{¥mathbf{¥theta}}¥[¥log{¥frac{p(x,¥mathbf{¥theta})}{p(x,¥mathbf{¥theta}’)}}¥]
  • 関数間の遠近関係
    • 関数は統計多様体上の点なので、関数の違いを多様体上の「距離」として測りたい
    • しかしながら、KL divergenceの定義は、「自身」とのdivergenceを0となるような定義式を使って、「他者」のdivergenceを正の値にしており、2つの関数のどちらを基準にするかで値が変わる。そのことを受けて「距離」ではなく"divergence"と呼ぶ
    • e-,m-flatsという2通りのパラメタ系の良さは次の通り
      • e-系で表しているとき、e-flatなので、ごちゃごちゃ考えずとも、e-系のパラメタを使って比較的簡単にdivergenceを計量することができる
      • m-系の場合も同じ
      • 今、3つの関数p,p',p''があったときに、
      • p,p'について、e-系表現がわかっていれば、両者のdivergenceが簡単に測れる
      • さらに、p',p''について、m-系表現がわかっていれば、両者のdivergenceも簡単に測れる
      • p,p''の間のdivergenceは、片やe-系、片やm-系なので、簡単に測れないけれど、ここで、e-系とm-系とが「直交」していることが役に立ち、両者が単純に足し合わせればよいことが示せる
  • 関数間の遠近を統計量にする
    • KL divergence(対数尤度の差)の2倍の値はカイ二乗分布に従うことが知られており、尤度比検定に持ち込める
  • Flatなsubmanifolds
    • e-,m-flatsのよいところは、それぞれの座標系で、関数空間全体がベクトル空間のようになっていることである
    • 今、例えばe-系で表された関数空間全体のうち、e-系の線形和で表せる部分空間(submanifold)を取り出したとする
    • これはe-系表現があるからe-flatだし、e-flatならばm-flatではあるけれど、m-系で線形和に表せるかどうかというとそうはならない
      • e-系とm-系とは、相互に「直交」関係にあり、e-のすべての要素のそれぞれをm-系の全要素が分担して表現する、というような関係になっている(なっていがち)なために、e-の方の部分をとっても、m-の方は「個々の要素の寄与分」は小さくなっても、mの要素数を減らすわけにはいかないことなどを念頭に置くとわかりやすい
  • うまく辻褄の合うe-系submanifoldシリーズとm-系submanifoldシリーズ
    • e-系とm-系とのパラメタには、対応関係があるので、それについて注意をしつつ、パラメタを除いていくと(パラメタの値を固定していくと)、入れ子になったe系submanifoldsシリーズとm-系submanifoldsシリーズができる。それは、直交性を利用して活用することができる

2016-05-19 90分で学ぶ医学・ライフサイエンスの統計学手法のいまどき

[][][][][][][][][]手法・目的の分類オーバービュー〜90分で学ぶ医学・ライフサイエンスの統計学手法のいまどき

  • Statistical Bioinformatics: For Biomedical and Life Science Researchers(の目次)
  • Quality Control of High-Throughput Data
  • Statistical Tests, Statistical Significance, Error Controlling
  • Classification/Clustering
    • Unsupervised Learning
    • Supervised Learning
  • Multidimensional/Highdimensional Data and Data visualization
  • Statistical Models and Inference
  • Random value-based approaches
    • Resampling
    • MCMC
  • Graph and Network
  • Experimental Designs

[][][]Multiple-Comparison Issue

  • When you test multiple times, you should not believe nominal p-values of individual tests. 検定を複数行った場合には、個々の検定のp値をそのまま使って解釈できない
p <- runif(10^5)
hist(p)
plot(sort(p),pch=20,cex=0.1)
alpha <- 0.05
abline(v=length(p)*alpha,col="red")
abline(h=alpha,col="blue")
  • Multiple tests (100 tests) 100個の検定
n <- 100
m <- 1000
ps <- matrix(runif(n*m),n,m)
min.ps <- apply(ps,2,min)
hist(min.ps)
plot(sort(min.ps),pch=20,cex=0.1)
alpha <- 0.05
abline(v=length(min.ps)*alpha,col="red")
abline(h=alpha,col="blue")
  • False Disovery Rate
    • q-value not p-value
    • Some tests should be "significant"; many genes should be significantly diffrently expressed between cancer and normal cells.
    • The smallest p-value should be small enough, but the 2nd, 3rd,... smallest could be not so small and they should be judged by different threshold.
n.null <- 80
n.alt <- 20
p <- runif(n.null)
p <- c(p,pchisq(rchisq(n.alt,df=1,ncp=10),df=1,lower.tail=FALSE))
p <- sort(p)
plot(p)
q <- p.adjust(p,"BH")
k <- max(which(q<0.05))
plot(p,pch=20)
abline(0,p[k]/k,col="red")

f:id:ryamada22:20160518140425j:image

plot(q,pch=20)

f:id:ryamada22:20160518140426j:image

    • Further reading : Large-scale inference; When you have many observations all together, use their overall distribution/behavior for judgement of each.

[]High Dimensionality

f:id:ryamada22:20160518124000p:image

  • Sparse
    • 10^5cells with 5 markers
      • 0.1 ¥times 10grid; one cubicle only has 1 cell in average. It is too sparse to estimate density in a regular way.
  • No center, no common individuals

f:id:ryamada22:20160518125810j:image

N <- 10^4
x <- runif(N)-0.5
y <- runif(N)-0.5

s <- which(x^2+y^2<0.5^2)

plot(x,y,pch=20,cex=1)
points(x[s],y[s],pch=20,cex=1,col="red")
n <- 1:10
r <- 0.5
vol.hypersphere <- pi^(n/2)/gamma(n/2+1) * r^n
plot(n,vol.hypersphere,pch=20,cex=3)

f:id:ryamada22:20160518150548j:image

[]small n large p

  • 100 samples x 25000 genes
  • You can predict perfectly when you are allowed to use explanatory variables as many as sample size.
n <- 10
m <- 10
p <- sample(0:1,n,replace=TRUE)
g <- matrix(sample(0:1,n*m,replace=TRUE),n,m)
p
g
lm.out <- lm(p~g)
p.est <- predict(lm.out)

plot(p,p.est)
print(cbind(p,p.est))
print(cbind(p,round(p.est)))

print(round(lm.out$coefficients))
> p
 [1] 1 1 0 0 0 1 1 0 0 0
> g
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    1    1    0    1    1    0    1    1     0
 [2,]    1    0    0    0    0    1    1    0    0     1
 [3,]    0    0    0    0    1    1    1    0    0     1
 [4,]    0    1    1    0    1    1    1    0    0     1
 [5,]    0    1    0    0    0    0    1    0    0     1
 [6,]    1    1    1    0    0    0    1    1    0     0
 [7,]    1    0    0    0    0    0    1    0    0     0
 [8,]    1    0    1    0    0    0    1    1    0     0
 [9,]    0    1    1    1    1    0    0    0    0     0
[10,]    1    0    0    0    1    1    1    0    0     0
> print(cbind(p,p.est))
   p         p.est
1  1  1.000000e+00
2  1  1.000000e+00
3  0 -3.333333e-01
4  0  8.881784e-16
5  0  3.333333e-01
6  1  6.666667e-01
7  1  6.666667e-01
8  0  3.333333e-01
9  0  1.110223e-16
10 0  3.333333e-01
> print(cbind(p,round(p.est)))
   p  
1  1 1
2  1 1
3  0 0
4  0 0
5  0 0
6  1 1
7  1 1
8  0 0
9  0 0
10 0 0
> 
> print(round(lm.out$coefficients))
(Intercept)          g1          g2          g3          g4          g5          g6          g7          g8 
          1           1           0           0          -1          -1           0          -1           0 
         g9         g10 
         NA          NA 
  • When you have p >> n explanatory variabls, you would find many variable sets, each of which perfectly predict dependent variable.

[]Noisy High-Throughput Biological Data

  • Noisy because
    • Biological phenomena are "noisy" ~ heterogeneity is the important feature of biology.
    • Experiments have many factors that add noise to data.
    • Highthroughput systems realize "highthoughput" by sacrificing preciseness somehow.
  • Noises are informative
    • Large-scale data have information on "shape" of noise and it can be used to quality-control/adjust the noisy data to purify "non-noisy signals".
  • 1000 samples, SNP chip for 50K SNP genotyping.
  • Success rate of individual samples.
N <- 1000
n1 <- 650
n2 <- 300

a <- rbeta(n1,50,2)
b <- rbeta(n2,100,20)
x <- rbeta(N-n1-n2,30,50)
hist(c(a,b,x),breaks = seq(from=0,to=1,length=40))

f:id:ryamada22:20160518145751j:image

[][]検定からベイズ推定へ Frequentist approach ~ tests and Bayesian approach

[][]MCMC

  • What is new?
  • This approach requires computers, that was not available in 20th century; that is why MCMC is being used for tasks that could not be solved by the methods before and also MCMC is being used for tasks that have been answered by the methods as well.
  • Link

[][][]時系列・空間とkrigingとか Time-series/Space and krigign etc.

  • Conventional analysis depends on independency among values
  • Time series/ spacial data have values that are next to each other are tightly associated and this feature should be absolutely considered.
  • At a glance

2016-05-18 Ricci flow と アルファ接続

[][][]Ricci flowとアルファ接続

  • ここ数日、統計多様体とその2つの平坦なパラメタの取り方について書いている
  • リーマン多様体として形を考えるときにRicci flowというのがあった
  • そこにもアルファ接続が出てくる
  • どういう関係なのか調べよう
  • こちらに資料