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

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

2017-06-18 Covariance Functions

[][][]Covariance Functions

  • 資料
  • Gaussian過程には分散共分散行列を使う
  • それをつくるには、2点に対して値を計算する関数が要る。それがCovariance functions
  • 2点の差ベクトルで決まると考える場合(x-x’)
  • 2点間距離で決まると考える場合(r=|x-x’|)
  • それ以外の距離っぽい定義で決まると考える場合(測地線距離、内積など)
  • などがある
  • rで決まるものが考えやすいが、rの関数はあまたある
  • 網羅的にそれを表すのが、Matern class of Covariance Functions
  • k_{Matern;¥nu}(r) = ¥frac{2^{1-r}}{¥Gamma(¥nu)} (¥frac{¥sqrt{2¥nu}r}{l})^¥nu ¥times K_¥nu (¥frac{¥sqrt{2¥nu}r}{l})
    • ただし、K_¥nu (x)= ¥frac{J_¥nu(x) ¥cos{(¥nu ¥pi)} - J_{-¥nu}(x)}{¥sin{(¥nu ¥pi)}} (第2種ベッセル関数)
    • ただし、J_¥nu (x) = ¥sum_{m=0} ^{¥infty} ¥frac{(-1)^m}{m! (m+¥nu +1)}(¥frac{x}{2})^{2m+¥nu} (第1種ベッセル関数)
    • これっておどろおどろしいのだけれど
      • ¥nu ¥to ¥inftyでSquare exponential covariance matrix(二乗距離に応じた指数関数的減少関数)
      • ¥nu = p + ¥frac{1}{2}を取らせると、指数関数とp次多項式で表される関数になる
        • k_{p+1/2}(r) = exp(-¥frac{¥sqrt{2¥nu}r}{l}) ¥frac{¥Gamma(p+1)}{¥Gamma(2p+1)} ¥sum_{i=0}^p ¥frac{(p+i)!}{i!(p-i)!} (¥frac{¥sqrt{8¥nu}r}{l})^{p-i}
      • ¥nuは微分可能な回数も定めていて、それは、Covariance functionに要請する滑らかさをコントロールする
      • 具体的には(それなりに滑らかにしたいならこの辺りがよい)
        • ¥nu = 3/2のとき(1+¥frac{¥sqrt{3} r}{l}) e^{-¥frac{¥sqrt{3} r}{l}
        • ¥nu = 5/2のとき(1+¥frac{¥sqrt{5} r}{l} + ¥frac{5 r^2}{3 l^2}) e^{-¥frac{¥sqrt{5} r}{l}
      • もっと粗くするなら¥nu = 1/2e^{-¥frac{r}{l}}。これは1次元なら酔歩のときの速度の変化モデルであるOrnstein-Uhlenbeck過程

2017-06-17 Gaussian Process Regression

[][]Gaussian Process Regression

  • こちらにGaussian Process Regressionをベタコードで解説してある
  • 手元のR環境だと、描図エラーが出るので、エラーを出す行だけコメントアウトしたものを以下に再掲する
  • xの座標ペアごとにどれくらいの関連を入れるかを分散共分散計算関数で指定し
  • 観測データと予測データとに関して、分散共分散行列を4分割したそれぞれの要素をGaussian Process的に作成し、それを線形代数的に解く、という過程を書いたコードになっている
  • ここが解いているところ
# These matrix calculations correspond to equation (2.19)
# in the book.
f.star.bar <- k.xsx%*%solve(k.xx)%*%f$y
cov.f.star <- k.xsxs - k.xsx%*%solve(k.xx)%*%k.xxs
# Demo of Gaussian process regression with R
# James Keirstead
# 5 April 2012

# Chapter 2 of Rasmussen and Williams's book `Gaussian Processes
# for Machine Learning' provides a detailed explanation of the
# math for Gaussian process regression.  It doesn't provide
# much in the way of code though.  This Gist is a brief demo
# of the basic elements of Gaussian process regression, as
# described on pages 13 to 16.


# Load in the required libraries for data manipulation
# and multivariate normal distribution
require(MASS)
require(plyr)
require(reshape2)
require(ggplot2)

# Set a seed for repeatable plots
set.seed(12345)

# Calculates the covariance matrix sigma using a
# simplified version of the squared exponential function.
#
# Although the nested loops are ugly, I've checked and it's about
# 30% faster than a solution using expand.grid() and apply()
#
# Parameters:
#	X1, X2 = vectors
# 	l = the scale length parameter
# Returns:
# 	a covariance matrix
calcSigma <- function(X1,X2,l=1) {
  Sigma <- matrix(rep(0, length(X1)*length(X2)), nrow=length(X1))
  for (i in 1:nrow(Sigma)) {
    for (j in 1:ncol(Sigma)) {
      Sigma[i,j] <- exp(-0.5*(abs(X1[i]-X2[j])/l)^2)
    }
  }
  return(Sigma)
}

# 1. Plot some sample functions from the Gaussian process
# as shown in Figure 2.2(a)

# Define the points at which we want to define the functions
x.star <- seq(-5,5,len=50)

# Calculate the covariance matrix
sigma <- calcSigma(x.star,x.star)

# Generate a number of functions from the process
n.samples <- 3
values <- matrix(rep(0,length(x.star)*n.samples), ncol=n.samples)
for (i in 1:n.samples) {
  # Each column represents a sample from a multivariate normal distribution
  # with zero mean and covariance sigma
  values[,i] <- mvrnorm(1, rep(0, length(x.star)), sigma)
}
values <- cbind(x=x.star,as.data.frame(values))
values <- melt(values,id="x")

# Plot the result
fig2a <- ggplot(values,aes(x=x,y=value)) +
  geom_rect(xmin=-Inf, xmax=Inf, ymin=-2, ymax=2, fill="grey80") +
  geom_line(aes(group=variable)) +
  theme_bw() +
  scale_y_continuous(lim=c(-2.5,2.5), name="output, f(x)") +
  xlab("input, x")


fig2a

# 2. Now let's assume that we have some known data points;
# this is the case of Figure 2.2(b). In the book, the notation 'f'
# is used for f$y below.  I've done this to make the ggplot code
# easier later on.
f <- data.frame(x=c(-4,-3,-1,0,2),
                y=c(-2,0,1,2,-1))

# Calculate the covariance matrices
# using the same x.star values as above
x <- f$x
k.xx <- calcSigma(x,x)
k.xxs <- calcSigma(x,x.star)
k.xsx <- calcSigma(x.star,x)
k.xsxs <- calcSigma(x.star,x.star)

# These matrix calculations correspond to equation (2.19)
# in the book.
f.star.bar <- k.xsx%*%solve(k.xx)%*%f$y
cov.f.star <- k.xsxs - k.xsx%*%solve(k.xx)%*%k.xxs

# This time we'll plot more samples.  We could of course
# simply plot a +/- 2 standard deviation confidence interval
# as in the book but I want to show the samples explicitly here.
n.samples <- 50
values <- matrix(rep(0,length(x.star)*n.samples), ncol=n.samples)
for (i in 1:n.samples) {
  values[,i] <- mvrnorm(1, f.star.bar, cov.f.star)
}
values <- cbind(x=x.star,as.data.frame(values))
values <- melt(values,id="x")

# Plot the results including the mean function
# and constraining data points

fig2a <- ggplot(values,aes(x=x,y=value)) +
  geom_rect(xmin=-Inf, xmax=Inf, ymin=-2, ymax=2, fill="grey80") +
  geom_line(aes(group=variable)) +
  theme_bw() +
  scale_y_continuous(lim=c(-2.5,2.5), name="output, f(x)") +
  xlab("input, x")




fig2b <- ggplot(values,aes(x=x,y=value)) +
  geom_line(aes(group=variable), colour="grey80") +
  #geom_line(data=NULL,aes(x=x.star,y=f.star.bar),colour="red", size=1) + 
  #geom_line(aes(x=x.star,y=f.star.bar),colour="red", size=1) + 
  geom_point(data=f,aes(x=x,y=y)) +
  theme_bw() +
  scale_y_continuous(lim=c(-3,3), name="output, f(x)") +
  xlab("input, x")

fig2b

# 3. Now assume that each of the observed data points have some
# normally-distributed noise.

# The standard deviation of the noise
sigma.n <- 0.1

# Recalculate the mean and covariance functions
f.bar.star <- k.xsx%*%solve(k.xx + sigma.n^2*diag(1, ncol(k.xx)))%*%f$y
cov.f.star <- k.xsxs - k.xsx%*%solve(k.xx + sigma.n^2*diag(1, ncol(k.xx)))%*%k.xxs

# Recalulate the sample functions
values <- matrix(rep(0,length(x.star)*n.samples), ncol=n.samples)
for (i in 1:n.samples) {
  values[,i] <- mvrnorm(1, f.bar.star, cov.f.star)
}
values <- cbind(x=x.star,as.data.frame(values))
values <- melt(values,id="x")

# Plot the result, including error bars on the observed points
gg <- ggplot(values, aes(x=x,y=value)) + 
  geom_line(aes(group=variable), colour="grey80") +
  #geom_line(data=NULL,aes(x=x.star,y=f.bar.star),colour="red", size=1) + 
  #geom_errorbar(data=f,aes(x=x,y=NULL,ymin=y-2*sigma.n, ymax=y+2*sigma.n), width=0.2) +
  geom_point(data=f,aes(x=x,y=y)) +
  theme_bw() +
  scale_y_continuous(lim=c(-3,3), name="output, f(x)") +
  xlab("input, x")

gg

2017-06-13 ぱらぱらめくる『Lecture notes on Bayesian Nonparametrics』

[][][]ぱらぱらめくる『Lecture notes on Bayesian Nonparametrics』

  • ぱらぱらめくれるか?
  • テキスト Lecture Notes on Bayesian Nonparametrics Peter Orbanz
  • 目次
    • Chapter 1 Terminology 用語・定義
    • Chapter 2 クラスタリングとディリクレ過程・チャイニーズ・レストラン過程
    • Chapter 3 Latent features/買い物記録とインディアン・ビュッフェ過程
    • Chapter 4 回帰とガウシアン過程
    • Chapter 5 モデルをピースとして組み立てる
    • Chapter 6 Exchangeability
    • Chapter 7 事後分布
    • Chapter 8 ランダム測度

[][][]Chapter 1 ノンパラベイズとは ぱらぱらめくる『Lecture notes on Bayesian Nonparametrics』

  • この章では、以下を確認する
    • パラメトリックとノンパラメトリック
    • (パラ・ノンパラの)ベイズ
  • パラメトリックとノンパラメトリック
    • 事象の空間がある(ベルヌーイなら{0,1}のこと、サイコロなら、{1,2,3,4,5,6}のこと、実数全体を事象空間とする確率分布には、正規分布ある、など)
    • ある事象空間に定め得る、ありとあらゆる確率質量/密度分布のことを考える。それらを「分布の集合」として考える
    • ありとあらゆるものの中には、「これこれの条件を満たす」というように定義できる分布だけをまとめると、それは部分集合
    • この部分集合を統計モデルと呼ぶ
    • 今、統計モデルの「これこれの条件を満たす」という部分を、パラメタで表すことにすれば、パラメタで表すことができる
    • このようにすることで、統計モデルは、パラメタが作る空間になり、個別の確率分布はこのパラメタが作る空間の点にな
    • ただし、このパラメタの数が、有限・一定な定義になっている場合もあれば、どこまででも増えてしまう場合もある。場合がある、というか、そのように定義できる
    • 有限個数のパラメタで定義づけした統計モデルをパラメトリックモデル、無限個のパラメタで定義づけした統計モデル(パラメタ数が無限個になりうる定義を持つ統計モデル)をノンパラメトリックモデルと言う
    • 有限個数パラメタモデルと無限個数パラメタモデルと呼ぶことにしてもよい
  • ベイズ
    • 統計モデルは有限個にしろ無限個にしろパラメタを持つ
    • パラメタの値について何かを考えたい(データからパラメタ値を決めたりすること)
    • ベイズとは、パラメタを確率変数として考えること
    • パラメタは何かしらの値を取りうるのだが、それを確率変数と捉え、事前分布・尤度・事後分布が定まる
    • 個々の分布は、パラメタの値のセットで決まるので、パラメタが確率変数であるから、個々の分布は、パラメタ空間の(分布に従う)ランダムな点であり、
    • 確率事象は、この「個々の分布」から発生するとみなしているから、もう一段、別のランダムなサンプリングを考える
  • ノンパラメトリック・ベイズ
    • ノンパラメトリック・ベイズでは、パラメタ空間が無限次元空間なので、無限次元空間に事前分布を定める必要がある
    • その方法がChapters 2,3,4で紹介され、それらを組み合わせる話が続く

[][][]Chapters 2,3,4 無限次元空間の分布を作る ぱらぱらめくる『Lecture notes on Bayesian Nonparametrics』

  • 無限次元分布発生方法の名前
    • ディリクレ過程、チャイニーズ・レストラン、Stick-breaking
    • インディアン・ビュッフェ過程
    • ガウス過程、Kriging、酔歩・ウィーナー過程、ガウシアンランダムフィールド
  • どういう推定課題なのかによる分類
  • Chapter 2 クラスタリング ディリクレ過程、チャイニーズ・レストラン、Stick-breaking
    • ノンパラ、無限個のパラメタ
    • クラスタ数が無限にあり、それらに、確率を持たせ、総和が1であるようにしたい
    • そんな無限クラスタへの確率の割付発生法の一つがStick-breaking。それをするときに、ベータ分布¥beta(1,¥alpha)を使うのがディリクレ過程
    • クラスタに割合が割り振られるが、そのクラスタが確率分布だとすると、その確率分布を決めるパラメタが必要。ディリクレ過程では、先ほどの¥alphaと、パラメタ発生させるための下となる分布が必要
    • クラスタリングを、標本の帰属クラスタIDの決定についてだけ考えるのであれば、「パラメタ値」は発生させなくてもよい。randam partitionと呼び、チャイニーズ・レストラン過程と言う
    • あるクラスタの割合をベータ分布で出し、次のクラスタの割合を、残り確率を1とみなしてその部分割合をベータ分布で指定して定めることを繰り返す
    • 別のやり方
      • あくまでも、クラス多数を無限大にして、クラスタの確率の総和を1にしたかっただけ。ディリクレ過程ではベータ分布を使った。それ以外の分布を使うこともある。特に、べき乗則に乗るような「割合の分散の大きな」クラスタ別頻度にするやり方もある
  • Chapter 3 特徴抽出 インディアン・ビュッフェ過程
    • ノンパラ、無限個のパラメタ
      • 特徴が無限個になりうる。個々の特徴に頻度を持たせると、パラメタ数は無限。さらに、個々の標本は複数の特徴を持ちうるので、特徴の組み合わせにも頻度を持たせることになり、パラメタ数は、特徴の個数 n ¥to ¥infty に対して¥lim_{n ¥to ¥infty}2^n-1自由度を持つ。
    • チャイニーズ・レストラン過程が
      • 既存のクラスタの1つに帰属させるか、新規のクラスタを1個作ってそれに帰属させるか、と言うプロセスなのに対して
      • 既存のクラスタのそれぞれに、帰属させるか否かを決め、さらに、新規のクラスタをポアソン分布に従って非負整数個生成し、それに追加サンプルを帰属させる、
      • という形式で特徴を増やしていく
  • Chapter 4 ガウス過程
    • ディリクレ過程、インディアンビュッフェ過程が無限個のクラスタ、無限個の特徴と言った「無限個の離散的な事象」を対象にしていたのに対して
    • ガウス過程は時空間と言った連続な台を対象にする
    • 連続な台なので、どうしても、「無限」に定める必要が出てきて、すなわち、無限次元〜ノンパラ
    • 関数を考える
    • 関数をd次元空間の点を1次元空間の点への対応付けとみなす
    • この関数について、d次元空間の点をn個取り出すことにすれば、それを1次元空間の(たかだか)n個の点に対応付ける、という関数を考えることができる
    • このような関数をランダムな存在とみなす
    • すると、n個の点の対応先である1次元空間の点はn次元空間のランダムな点になる
    • そのn次元空間のランダムな点を確率密度分布とみることもできる
    • ここで、このランダムな点をn次元正規分布と仮定しよう、というのがガウス過程
    • これは、d次元空間の点に対応する1次元空間の点の座標には、期待値と分散があるということだし、d次元空間の点のペアに対応する1次元空間の点の座標のペアには共分散があるということ
    • n次元正規分布なので、期待値(中心座標)と分散共分散行列が定まる
    • 近い点同士に対応づく1次元空間座標が近いべきである、という気持ちは共分散の大きさで表され、それは関数で言えば滑らかさに対応する
    • ガウス過程は、関数の滑らかさを分散共分散行列で表現しようというモデルであるとも言える
      • より詳しく言うと、分散共分散行列は stationarity, isotropy, smoothness and periodicityを定める(らしい)Wiki
    • 分散共分散行列の分布に事前分布を入れることで、関数推定が線型代数処理を通じて得られる
    • ガウス過程が仮定しているのは、期待値ベクトルと分散共分散行列。そこに何がしかの観察がなされると、未観察に関して、事後分布が得られるが、その事後分布は共分散によって観察値に引きずられた分布になる。その事後分布は、仮定している分散共分散行列と、実際に観察した値が指し示す分散共分散(その中にランダムエラーも入る)との線型計算
    • 酔歩、ガウシアン・ランダムフィールドが含まれる
    • 分散共分散行列にいろいろあって、それらが、stationarity, isotropy, smoothness and periodicityの定義と対応づくWikiにそれが書いてあるようだ
  • Chapter 5
    • 組み合わせる
    • (無限にある)離散項に重みをつけて、その重み付き和を取ることが混合モデルだった
    • 連続に拡張すればそれは重み関数付きの積分
    • それは関数の掛け算
    • 関数の掛け算は『無限項』に相当するが、関数がパラメタライズされていれば、そのパラメタのハンドリングに簡略化される
  • Chapter 6 Exchangeability
    • 標本がiidなとき、その標本の順序が変えられる
    • 無限次元・ノンパラになると、Exchangeabilityとiidとが対応づくわけではなくなるらしい
    • 条件付き確率で定義されるらしい

2017-06-10 Rstan 基礎

[][][]Stanサンプリング

---
title: "StanDistribution"
author: "ryamada"
date: "2017年6月10日"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rstan)
```


# Stanの各種分布の練習

```{r}
.stan_code1 = '
data {
  int n_obs;
  vector [n_obs] x;
}

parameters {
  real mu;
  real<lower=0> sigma;
}

model {
  x ~ normal(mu, sigma);
}'
```
```{r,echo=FALSE}
#Compiling
MODEL1 = stan_model(model_code = .stan_code1)
```

```{r,echo=FALSE}
.stan_code2 = '
data {
  int n_obs;
  vector [n_obs] x;
}

parameters {
  real mu;
  real<lower=0> sigma;
}

model {
  x ~ cauchy(mu, sigma);
}'
#Compiling
MODEL2 = stan_model(model_code = .stan_code2)
```
```{r,echo=FALSE}
.stan_code3 = '
data {
  int n_obs;
  vector [n_obs] x;
}

parameters {
  real a;
  real b;
}

model {
  x ~ gamma(a,b);
}'
#Compiling
MODEL3 = stan_model(model_code = .stan_code3)
```
```{r,echo=FALSE}
.stan_code4 = '
data {
  int n_obs;
  int x[n_obs];
}

parameters {
  real <lower=0> a;
}

model {
  x ~ poisson(a);
}'
#Compiling
MODEL4 = stan_model(model_code = .stan_code4)
```

```{r,echo=FALSE}
.stan_code5 = '
data {
  int n_obs;
  int x[n_obs];
  int ns[n_obs];
}

parameters {
  real <lower=0> p;
}

model {
  x ~ binomial(ns, p);
}'
#Compiling
MODEL5 = stan_model(model_code = .stan_code5)
```
```{r, echo=FALSE}
.stan_code6 = '
data {
  int n_obs;
  int K;
  int<lower=1,upper=K> Z[n_obs];  # data: category index
  vector<lower=0>[K] alpha;
}

parameters {
  simplex[K] theta; 
}

model {
  for (i in 1:n_obs)
   Z[i] ~ categorical(theta);
  theta ~ dirichlet(alpha); 
}'
#Compiling
MODEL6 = stan_model(model_code = .stan_code6)

```
```{r,echo=FALSE}
.data1 = list(x=rnorm(10000, 10, 3))
.data1$n_obs = length(.data1$x)
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit1 = sampling(MODEL1, data=.data1)


EXT1 = extract(.fit1)
```
```{r,echo=FALSE}
.data2 = list(x=rcauchy(10000, 10, 3))
.data2$n_obs = length(.data2$x)
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit2 = sampling(MODEL2, data=.data2)


EXT2 = extract(.fit2)
```
```{r,echo=FALSE}
.data3 = list(x=rgamma(10000, 2,5))
.data3$n_obs = length(.data3$x)
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit3 = sampling(MODEL3, data=.data3)


EXT3 = extract(.fit3)
```
```{r,echo=FALSE}
.data4 = list(x=rpois(10000, 2.5))
.data4$n_obs = length(.data4$x)
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit4 = sampling(MODEL4, data=.data4)


EXT4 = extract(.fit4)
```

```{r,echo=FALSE}
.data5 = list(ns=rpois(10000, 2.5))
p = 0.3
.data5$x = rbinom(length(.data5$ns),.data5$ns,p)
.data5$n_obs = length(.data5$x)

#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit5 = sampling(MODEL5, data=.data5)


EXT5 = extract(.fit5)
```
```{r,echo=FALSE}
K <- 5
library(MCMCpack)
r <- rdirichlet(1,rep(1,K))
.data6 = list(Z=sample(1:K,10000,replace=TRUE,r))
.data6$K = K
.data6$n_obs = length(.data6$Z)
.data6$alpha = rep(1,K)

#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit6 = sampling(MODEL6, data=.data6)


EXT6 = extract(.fit6)
```
```{r}
print(.fit1)
summary(.fit1)
plot(.fit1)
pairs(.fit1)
rstan::traceplot(.fit1)
rstan::stan_trace(.fit1)
rstan::stan_hist(.fit1)
rstan::stan_dens(.fit1)
```
```{r}
print(.fit2)
summary(.fit2)
plot(.fit2)
pairs(.fit2)
rstan::traceplot(.fit2)
rstan::stan_trace(.fit2)
rstan::stan_hist(.fit2)
rstan::stan_dens(.fit2)
```
```{r}
print(.fit3)
summary(.fit3)
plot(.fit3)
pairs(.fit3)
rstan::traceplot(.fit3)
rstan::stan_trace(.fit3)
rstan::stan_hist(.fit3)
rstan::stan_dens(.fit3)
```
```{r}
print(.fit4)
summary(.fit4)
plot(.fit4)
pairs(.fit4)
rstan::traceplot(.fit4)
rstan::stan_trace(.fit4)
rstan::stan_hist(.fit4)
rstan::stan_dens(.fit4)
```
```{r}
print(.fit5)
summary(.fit5)
plot(.fit5)
pairs(.fit5)
rstan::traceplot(.fit5)
rstan::stan_trace(.fit5)
rstan::stan_hist(.fit5)
rstan::stan_dens(.fit5)
```

```{r}
print(.fit6)
summary(.fit6)
plot(.fit6)
pairs(.fit6)
rstan::traceplot(.fit6)
rstan::stan_trace(.fit6)
rstan::stan_hist(.fit6)
rstan::stan_dens(.fit6)
```

[][]データ型

  • [:title=参考]
  • それぞれの型にRリストでデータを渡す
  • データを渡すだけで、Stan処理としては、ただの1次元正規分布のパラメタ推定
  • stanでは、ベクトルを要素とする行列、とか、行列を要素とするベクトルとかがあるが、それは、Rのアレイで渡せばよいようだ
---
title: "StanDataType"
author: "ryamada"
date: "2017年6月10日"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rstan)
```

# データタイプ

```{r}
.stan_code1 = '
data {
  # スカラー
  real r;
  int n;
  int m;
  # ベクトル、行列
  vector[n] vc;
  row_vector[m] vr;
  matrix[n,m] M;
  # アレイ
  real x[n];
  int y[n];
  # 要素がベクトル・行列である配列
  vector[m] v[n];
  matrix[n,m] K[n];
  # 次元一般化アレイ
  real A[n,m];
  real B[n,n,m];
  # 要素が行列である行列
  matrix[n,m] C[n,m];

}

parameters {
  real mu;
  real<lower=0> sigma;
}

model {
  x ~ normal(mu, sigma);
}'
```
```{r,echo=FALSE}
#Compiling
MODEL1 = stan_model(model_code = .stan_code1)
```

```{r,echo=FALSE}
n <- 10
m <- 5
.data1 = list(x=rnorm(n, 10, 3))
.data1$n = n
.data1$m = m
  # スカラー
r = 0.4
.data1$r = r

.data1$vc = rnorm(n)
.data1$vr = rnorm(m)
.data1$M = matrix(rnorm(n*m),ncol=m);
.data1$y = rpois(n,r)
  # 要素がベクトル・行列である配列
.data1$v = matrix(rnorm(n*m),ncol=m)
.data1$K = array(rnorm(n*n*m),c(n,n,m))
  # 次元一般化アレイ
.data1$A = matrix(rnorm(n*m),ncol=m);
.data1$B= array(rnorm(n*n*m),c(n,n,m))
  # 要素が行列である行列
.data1$C = array(rnorm(n*m*n*m),c(n,m,n,m))
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit1 = sampling(MODEL1, data=.data1)


EXT1 = extract(.fit1)
```

2017-06-04 『法数学入門』を書く

[][][]『法数学入門』を書く

  • まずは章立てを検討する
---
title: "法数学入門"
author: "ryamada"
date: "2017年6月4日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
library(MCMCpack)
library(rgl)
```

# はじめに

# 基本的数学
## 関数、指数関数・対数関数
## 微分・積分
## 順列・組み合わせ
## 法数学ではこう使う

# 確率と尤度
## 確率関数と尤度関数
## 二項分布とベータ分布
## 多項分布とディリクレ分布
## 尤度比
## 法数学ではこう使う
### 二者択一の仮説の比較
### 多アレルの遺伝的マーカー

# 独立と非独立
## 独立の仮定と確率・尤度
## 非独立の仮定と確率・尤度
## 法数学ではこう使う
### 非血縁関係と血縁関係
### 連鎖と連鎖不平衡

# たくさんあること
## 多重検定補正
## 法数学ではこう使う
### 誕生日のパラドクス
### 70億人

# ベイズ推定
## 最尤推定とモデル
## 事後分布推定と点推定と区間推定
## 法数学ではこう使う

# グラフィカルモデル、ベイジアンネットワーク
## 機械学習
## ベイジアンネットワーク
## 法数学ではこう使う

# コンピュータアプリケーションとコンピュータ言語
## 統計・シミュレーションとR
## 機械学習とPython