2010-02-01から1ヶ月間の記事一覧

平均0、分散1の正規分布は 自由度kのカイ分布(カイ自乗ではなく)は k次元球の表面積は これらを使うと、自由度kのカイ分布は これは、が表すように、原点から、遠ざかると確率が小さくなる分布であって、その小さくなり方が、正規分布と同じように、が一…

一般正規分布・一般誤差分布・指数べき分布

変数の置き方を変えると とか とか書ける。 のときは指数分布を正負で対象にしてやったもの のときは正規分布 これは、であることによる。

分布をすこしずつ変える

今日(2月25日)は国公立大学の入学試験。春のような暖かさの中、緊張した顔つきの受験生が見えます。緊張することはいいことです。みんな、がんばって欲しいです。 わからないことがたくさんあって、それがわかることは素敵ですから、是非、大学へ! さて…

矩形分布から指数分布・正規分布を経て一様分布へ(誤差分布)

2dhist()

n1 <- 500 n2 <- 300 n3 <- 200 x <- c(rnorm(n1, 0, 0.5), rnorm(n2, 5, 1), rnorm(n3, 8, 2)) y <- c(rnorm(n1, 0, 2), rnorm(n2, 3, 2), rnorm(n3, -3, 1)) library(gregmisc) plot(x, y) h2d <- hist2d(x, y, show = FALSE, same.scale = TRUE, nbins = …

Rで2次元ヒストグラム

MDS

サンプルの多次元データが与えられているときに、特異値分解や、分散共分散行列の固有値分解をする話しが昨日まで。 MDSっていうのもある。 こちらは、サンプルに関してペアワイズの距離が与えられているとき、サンプル間の内積行列を再構成してやって、それ…

MDS 多次元尺度構成法と固有値分解

中心化後特異値分解と固有値分解

参考こちら が特異値分解。 変形して 。 これを解いて、S,V->Uが得られる。 今、Xを中心化すると はXの分散共分散行列に比例した値になるので、中心かした特異値分解と固有値分解は、同じようなもの。 #構造化集団をシミュレート Nm<-1000 #マーカー数 Npop<…

ジェノタイプデータのPCA補正

昨日の続き PCAにより、GWASのジェノタイプデータでいくつかの軸情報で個人に「位置情報」が与えられる 個人の位置情報に応じて、個人のフェノタイプ(ケース・コントロール)とSNPのジェノタイプの値を補正する 補正したフェノタイプと補正したジェノ…

RでSNPデータのEigenstrat検定補正

ジェノタイプデータのPCAその2

昨日の続きでは、正方行列を作らずに、非正方行列のままsvd()をかけるとどうなるかもやってみます。 同じ構造を表す固有値と固有ベクトルが取れました。 #構造化集団をシミュレート Nm<-1000 #マーカー数 Npop<-4 #亜集団数 Ns<-c(100,150,200,250) #集団別…

RでSNPデータのPCAプロットその2 svd()

ジェノタイプデータのPCA eigen()

集団構造化があるときに、PCAして、プロットすることがある。 その情報を使って、形質マッピング検定に用いる前座のようなもの。 この論文がEigenstratのそれですが。 ここでやっている、PCA部分をRでなぞってみます #構造化集団をシミュレート Nm<-1000 #マ…

RでSNPデータのPCAプロット eigen()

Rのimage()関数を使ってペアワイズLDプロットを簡単に描く方法を以前書いた(こちら)。 今日は、ハプロタイプの01表記ファイル"hoge.txt"から描くことにする。 横がマーカー、縦に染色体。こんな感じ。 0 0 0 0 1 0 1 1 0 1 0 1 1 0 0 1 1 1 0 0 1 0 0 0 1 …

いろいろ数え上げる

数え上げ # 順列 permN<-function(N=10,k=3){ return(exp(lgamma(N+1)-lgamma((N-k)+1))) } 組み合わせ combN<-function(N=10,k=3){ return( exp(lgamma(N+1)-lgamma((N-k)+1)-lgamma(k+1)) ) } 重複順列 repPermN<-function(N=10,k=3){ return(N^k) } 重複…

RでLDプロット その2

ROCカーブのAUCの信頼区間に関するメモ

ROC,AUCに関する概論はこちら。 Rには、DiagnosisMed というパッケージがあって、それにROC()という関数がある。 出力はこんな感じ Sample size: 170 Sample prevalence: 0.4118 Population prevalence: 0.4118 - same as sample prevalence if not informed…

アレルの分岐図

家系図は、個人の遺伝的伝達を図示したもの 染色体上のアレルも伝達される。 いくつも木ができる 染色体数を固定する。すべての染色体がペアを作って、そのペアに2本の新たな染色体を作らせる経過で、木がどうなるかを描かせてみる ローカスの移動に伴って…

Rでコアレセント シミュレーション

記号

家系図記号というのがある。だれがどこで統一指標を出しているのかと気になった。 日本人類遺伝学会のサイトには(すぐには)見当たらず、信州大学のサイトには、「アメリカ人類遺伝学会」が参照先として記載されていた(こんな感じ)。 アメリカ人類遺伝学会の…

クラスタリングの場合の数 その2

Rにphangornというパッケージ(CRAN記事はこちら)があって、その中のallTrees()という関数があって、それが数え上げてくれる。ただし、すべての木を作る関数なので、葉の数は10個まで。 library(phangorn) allTrees(5) 15 phylogenetic trees trees <- all…

クラスタリングの場合の数

今、N個のサンプルがあって、これを2分岐木のクラスタに纏め上げたいとする 何通りの木の形状(トポロジー)があるんだろう? 漸化式で考える がその数とする である のときを考える まではわかっているとして、Nの場合を知りたいものとする はに分けることが…

2分岐木のトポロジーを数え上げる?

evd 極値分布パッケージ

繰り返し観察の極端な値(最大1時間降水量とか。多重検定の最大統計量とかも・・・)の確率密度分布は極値分布と呼ばれる分布になることが知られている その形は、一般に、3つのパラメタ(location:,scale:,shape:)で規定される、以下の式で表される この一般式…

Rの極値分布パッケージでマルチプルテスティング

ハッセ図

要素数kの集合のべき集合は空集合を含めての部分集合からなり、それらの包含関係は、半順序になる。 半順序の関係をグラフに表すときにハッセ図を用いることがあるが、それを、Rで描いてみよう。 powerSetHasse<-function(k=4){ lev<-k+1 numelem<-rep(0,le…

Rでハッセ図