2012-07-26

Rを使ってみる

MIKU R

Rの操作を覚えてみる、の会
ベクトルを扱った
ベクトルの要素同士の加減乗除ができたら…
内積が出せる
内積が出せれば、ベクトルの長さも出せる
内積が出せれば、三角形の面積が出せる
- 三角形の面積の公式はこんなにある
三角形は台形の特殊形とみなす話がこちらにある
- ちょっと書いてみたけれど、高次元台形の位置座標の取扱いとか、そこに『平行線』の端点である頂点の対応とかを入れ込んだオブジェクトを作るのが面倒くさいことに気付いた
- 気づくまでに適当に書き散らしたソースを貼っておく
  - うまく回らないままです
もう一つ考えた
台形の下底を固定して、上底の端点を決める。台形というのは、上底の端点からの辺が上底と水平な特殊な場合で、それ以外は、一般の四角形になる、とか、上底をベクトル扱いする(向きを考慮する)と、面積は、ただしは上底と下底の長さ、というより、向き考慮の平均とかになりそうだった
- そのときの面積が幾何的にどうなるのか、とかも興味深かった
- じゃあ、というので、端点を固定して、「上底になるかもしれない辺」の長さも固定して、上底になるかもしれない辺の方向をぐるりと１周させるっていうのは、どういうことなのか、とかもやってみた(後半のソース)

# n次元台形

n <- 3

# n-2次元空間の点がn個

X <- matrix(rnorm((n-1)*n),ncol = n)

# n個のn-1次元空間の点のそれぞれについて、線分をとり、
# その線分はすべて並行

Y <- matrix(rnorm(n*2),ncol = n)

# n次元台形の頂点座標を列ベクトルとする行列
Z <- cbind(rbind(X,Y[1,]),rbind(X,Y[2,]))

Z


# 別の方法でn次元台形

n <- 3

# 台形の平行辺の方向単位ベクトル
v <- runif(n)
v <- v/sqrt(sum(v^2))

# n個のn次元空間の点

X <- matrix(rnorm(n^2),ncol = n)

# それぞれの点からの平行辺のもう片方の端への方向付きベクトル長

R <- rnorm(n)

Y <- X + v %*% t(R)

calc.theta <- function(v){
	n <- length(v)
	theta <- rep(0,n)
	tmp <- v
	for(i in n:1){
		theta[i] <- asin(tmp[i])
		tmp[1:i] <- tmp[1:i]/cos(theta[i])
	}
	theta
}
t1 <- pi/3
t2 <- pi/4
v2 <- c(cos(t1)*cos(t2),sin(t1)*cos(t2),sin(t2))

calc.theta(v2)

make.mat <- function(thetas){
	n <- length(thetas)
	ret <- diag(n)
	for(i in 2:n){
		tmp <- diag(n)
		tmp[1,1] <- tmp[i,i] <-cos(thetas[i])
		tmp[1,i] <- -sin(thetas[i])
		tmp[i,1] <- -tmp[1,i]
		ret <- ret %*% tmp
	}
	ret
}

make.rot2v <- function(v){
	thetas <- calc.theta(v)
	m <- make.mat(thetas)
	return(list(m = m,thetas = thetas))
}


my.hv <- function(X,Y){
	n <- length(X[,1])
	if(n == 0){
		return(1)
	}
	P <- X-Y
	v <- apply(P,1,sum)
	v <- v/sqrt(sum(v^2))
	m <- make.rot2v(v)$m
	print(m)
	X.2 <- t(m) %*% X
	Y.2 <- t(m) %*% Y
	
	L <- X.2[1,]-Y.2[1,]
	ret <- sum(L)/n * my.hv(X.2[2:n,])
	return(ret)
}

my.hv(X,Y)
# n次元台形

n <- 3

X <- list()

X[[1]] <- matrix(0 ,nrow = 1, ncol =1)
for(i in 1:n){
	tmp1 <- rbind(X[[i]],rnorm(length(X[[i]][,1])))
	tmp2 <- rbind(X[[i]],rnorm(length(X[[i]][,1])))
	X[[i+1]] <- cbind(tmp1,tmp2)
}

X

triangle.area <- function(x,y){
	1/2 * sqrt(sum(x^2)*sum(y^2) - sum(x*y)^2)
}

L1 <- 1

L2 <- 1

x2 <- 0
y2 <- 1

theta <- seq(from = 0, to = 1, length = 10000) * 2 *pi
theta <- theta[-length(theta)]
x0 <- x2 -y2/tan(theta)

x3 <- x2 + L2 * cos(theta)
y3 <- y2 + L2 * sin(theta)

tri1 <- tri2 <- rep(0,length(theta))
for(i in 1:length(theta)){
	tri1[i] <- triangle.area(c(0,0)-c(x0[i],0),c(x2-x0[i],y2))
	tri2[i] <- triangle.area(c(L1,0)-c(x0[i],0),c(x3[i]-x0[i],y3[i]))
}

tri2_1 <- tri2-tri1

plot(theta,(tri2_1))

abline(h = 1/2 * y2*(L1+L2))
abline(h = 1/2 * y2*(L1-L2))