- Rmdファイルが以下にあります
- Rmdファイルからhtmlファイルを作る、epubファイルを作るには、こちらを参照
---
title: "一様分布乱数 Random values in uniform distribution"
author: "ryamada"
date: "Thursday, January 29, 2015"
output: html_document
---
p=1次元の単位区間[0,1]の一様分布の生成。
Generation of random values in p(=1)-dimentional unit segment [0,1].
```{r}
library(knitr) # TO embed 3d plots in the html documents
n <- 1000 # No. random varlues
R1 <- runif(n)
hist(R1)
```
p次元単位区間はp=2ならば正方形、p=3ならば立方体、一般にp次元超立方体。
p-dimensional unit "segments" are square for p=2, cube for p=3 and in general p-hypercubes.
```{r}
p <- 3
Rp <- matrix(runif(n*p),ncol=p)
#install.packages("rgl") # package to make 3D plot
library(rgl)
```
```{r setup}
knit_hooks$set(rgl = hook_rgl)
```
```{r,rgl=TRUE}
plot3d(Rp)
```
p=1次元の標準正規乱数を発生する。
Generate random values in p=1-dimensional standard normal distribution.
```{r}
Rn <- rnorm(n)
hist(Rn)
```
p次元正規乱数を発生する。
Generate random values in p-dimensional normal distribution.
```{r}
p <- 2
Rn2 <- matrix(rnorm(n*p),ncol=p)
plot(Rn2,asp=1) # "asp=1" to make scale of x and y axes the same
```
```{r, rgl=TRUE}
p <- 3
Rnp <- matrix(rnorm(n*p),ncol=p)
plot3d(Rnp)
```
p次元正規分布は中心からの距離に依存し、方向によらない。全方向について平等な分布。
p-normal distributions only depends on distance from the center not on directions from the center, meaning equidirectional.
単位円周上の一様乱数を発生する。
Generate random values in uniform distribution in the unit circle.
中心からの距離が1で角度が[0,2pi]について一様。
Distance from the center is 1 and angle is in uniform distribution in [1,2pi].
```{r}
r <- rep(1,n)
theta <- runif(n)*2*pi
Rc <- cbind(r*cos(theta),r*sin(theta))
plot(Rc,asp=1)
```
３次元以上だと困る。
This is not good for p >=3.
正規分布の方向に関する平等性を使う。
Use equidirectional feature of normal distribution.
```{r}
p <- 2
Rc2 <- matrix(rnorm(n*p),ncol=p)
plot(Rc2,asp=1) # direction-independent
```
これを使って円周上の乱点を作る。
Using this random points, make random points in the unit curcle.
```{r}
# calculate distance from the center
d <- sqrt(apply(Rc2^2,1,sum))
# Make distance of points being 1
Rc2. <- Rc2/d
plot(Rc2.,asp=1)
```
３次元以上でも大丈夫。
This is good for p>=3.
```{r, rgl=TRUE}
p <- 3
Rcp <- matrix(rnorm(n*p),ncol=p)
d <- sqrt(apply(Rcp^2,1,sum))
Rcp. <- Rcp/d
plot3d(Rcp.)
```
円周だけでなく、内側も含めた「円板」の一様分布。
The uniform distribution not only in the circle edge, but including the inside of the circle disc.
円板を含む正方形の一様分布から、円板部分を取り出す。
Uniform distribution of unit square first, then take out the disc area.
```{r}
p <- 2
Rp <- matrix(runif(n*p,-1,1),ncol=p)
plot(Rp,asp=1)
```
```{r}
d <- sqrt(apply(Rp^2,1,sum))
s <- which(d<=1)
Rs <- Rp[s,]
plot(Rs,asp=1)
```
n個の乱点を発生するには、多めに乱点を発生させて、円板上の点をn個確保する必要がある。
When you need n random points in the disc, generate more than n points in the square first, then select n points in the disc.
```{r}
N <- n*10 # Larger number to generate first
p <- 2
Rp <- matrix(runif(N*p,-1,1),ncol=p)
plot(Rp,asp=1)
d <- sqrt(apply(Rp^2,1,sum))
s <- which(d<=1)
length(s)
```
```{r}
s. <- sample(s,n) # select n points
length(s.)
Rs <- Rp[s,]
plot(Rs,asp=1)
```
p=2のときは円板、p=3のときは、中身の詰まった単位球。
p=2 unit sphere disc is a circle with inside, p=3 unit sphere disc is a sphere with inside.
```{r}
N <- n*10 # Larger number to generate first
p <- 3
Rp <- matrix(runif(N*p,-1,1),ncol=p)
plot(Rp,asp=1)
d <- sqrt(apply(Rp^2,1,sum))
s <- which(d<=1)
length(s)
```
```{r, rgl=TRUE}
s. <- sample(s,n) # select n points
length(s.)
Rs <- Rp[s,]
plot3d(Rs,asp=1)
```
pを増やして、単位p次元円板内の点の数を数える。
Make p bigger and count number of points inside p-dimensional unit sphere disc.
```{r}
ps <- 2:10
cnt <- rep(0,length(ps))
N <- 10^4
for(i in 1:length(ps)){
p <- ps[i]
Rp <- matrix(runif(N*p,-1,1),ncol=p)
d <- sqrt(apply(Rp^2,1,sum))
s <- which(d<=1)
cnt[i] <- length(s)
}
cnt
plot(ps,cnt)
```
正規分布は、方向に関して平等。
Normal distribution is equi-directional.
それを用いて、球面一様分布を発生させた。
Using it, uniform distribution in unit sphere surface wa s generated.
単位球面上の点は中心からの距離がすべて1.
The distance of the ooints on the sphere from the center is 1.
中心からの距離が[0,1]の点がほしい。
We want to make points with distance in [0,1].
p=2次元円板の場合、中心からの距離がxの点は長さ$2\pi x$の周にあるから、中心からの距離xに比例して、点の数が増えるはずである。
今、中心からの距離xをその値に比例して増えるように与えたい。
Consider p=2-dimensional circle disc. Points whose distance from the center is x is on a circle whose length is $2/pi x$. Therefore number of points should increase proportinal to x.
$pdf(x) = \frac{1}{2}x$ をその確率密度分布である。
その積分が累積分布であり、それは$cdf(x) = x^{2}$である。
確率密度分布とその累積分布を描いておく。
Its probability density function(pdf) and cumulative distribution function(cdf) are $pdf(x) = \frac{1}{2}x$ and $cdf(x) = x^{2}$, respectively. We can draw them.
```{r}
p <- 2
x <- seq(from=0,to=1,length=1001)
pdf <- 1/p* x^{p-1}
cdf <- x^{p}
matplot(x,cbind(pdf,cdf),type="l")
```
累積分布の逆関数 $x = q^{\frac{1}{2}}$を使えば、短い方から割合(クオンタイル）がq であるような距離xが求められる。
Using the inverse function of cdf(x), $x = q^{\frac{1}{2}}$, we can calculate distance from the center for arbitrary quantile values.
まず、p=2次元円周一様分布を描く。
Draw uniform random points in p=2-unit circle, first.
```{r}
p <- 2
Rcp <- matrix(rnorm(n*p),ncol=p)
d <- sqrt(apply(Rcp^2,1,sum))
Rcp. <- Rcp/d
plot(Rcp.)
```
ついで、この点の半径を[0,1]の範囲で、半径に比例した確率分布でランダムに与える。
Next, randomly assign radius to each point, where the radius is in distance-proportinal density function.
```{r}
q <- runif(n) # uniform random values as quantile values
Q <- q^(1/(p)) # Inverse function of cdf
Rcp.. <- Rcp. * Q
plot(Rcp..)
```
p=3
中心からの距離xに関して、$x^{p-1}$に比例して点の数は増える。
Number of points increase proportinal to $x^{p-1}$, where x is distance from the center.
pdf and cdf.
```{r}
p <- 3
x <- seq(from=0,to=1,length=1001)
pdf <- 1/p* x^{p-1}
cdf <- x^{p}
matplot(x,cbind(pdf,cdf),type="l")
```
単位球面一様分布。
Uniform distribution on sphere surface.
```{r, rgl=TRUE}
p <- 3
Rcp <- matrix(rnorm(n*p),ncol=p)
d <- sqrt(apply(Rcp^2,1,sum))
Rcp. <- Rcp/d
plot3d(Rcp.)
```
中心からの距離を与える。
Give radius to points.
```{r, rgl=TRUE}
q <- runif(n) # uniform random values as quantile values
Q <- q^(1/(p)) # Inverse function of cdf
Rcp.. <- Rcp. * Q
plot3d(Rcp..)
```
make.kame.hyoushi <- function(ttl="",n.kame=30,kame.col="blue",sub=""){
t <- seq(from=0,to=1,length=1000)*2*pi
kora <- cbind(cos(t),sin(t))
s <- 1/4
teashi <- kora*s
thetas <- c(0,1,2,4,5)/6*2*pi
teashi.ctr <- cbind(cos(thetas),sin(thetas))
teashi.coords <- matrix(0,nrow=0,ncol=2)
for(i in 1:length(thetas)){
teashi.coords <- rbind(teashi.coords,cbind(teashi.ctr[i,1]+teashi[,1],teashi.ctr[i,2]+teashi[,2]))
}
teashi.out <- teashi.coords[which(apply(teashi.coords^2,1,sum)>1),]
shippo.st <- c(cos(pi),sin(pi))
shippo.s <- 0.1
shippo.l <- 0.4
shippo.x <- seq(from=0,to=-shippo.l,length=100)
shippo.y <- shippo.s * sin(shippo.x*2*pi/shippo.l)
shippo.coords <- cbind(shippo.st[1] + shippo.x,shippo.st[2]+shippo.y)
kame <- rbind(kora,teashi.out,shippo.coords)
b <- 0.9
kame2 <- kame
kame2[,2] <- kame2[,2]*b
X <- matrix(0,nrow=0,ncol=2)
for(i in 1:n.kame){
theta <- runif(1)*2*pi
R <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),2,2)
new.kame <- t(R %*% t(kame2))
mv <- runif(2)*n.kame/2
new.kame <- cbind(new.kame[,1]+mv[1],new.kame[,2]+mv[2])
X <- rbind(X,new.kame)
}
jpeg(filename = paste(ttl,".jpeg"), width = 5000, height = 5000,
quality = 100, pointsize = 180,bg = "white", res = NA,
restoreConsole = TRUE)
plot(X,col=grey(1),axes=FALSE,xlab="",ylab="",main=paste(ttl,"\n",sub))
points(X,cex=0.5,pch=20,xlim=range(X),ylim=range(X),col=kame.col)
dev.off()
}
n.kame <- 8
kame.col <- "orange"
ttl <- "一様分布：超立方体・超球面・超円板"
subttl <- ""
make.kame.hyoushi(ttl,n.kame,kame.col,subttl)
