■-外資系企業勤務の奮闘記-■

R

R の日本語文章 (pdf 版)
http://buran.u-gakugei.ac.jp/~mori/LEARN/R/

R on Windows
http://plaza.umin.ac.jp/~takeshou/R/

R - Tips　- 統計解析ソフト R の備忘録頁 -
http://cse.naro.affrc.go.jp/takezawa/r-tips/r2.html

とりあえず，第1弾として，シミュレーションを行うにあたり，最低限の基本プログラムです．
#－－－－－－COINの乱数(1回発生させる場合)－－－－－－
coin<-function()
{
x<-runif(1)
if (x<=1/2) men<-1 #表
else men<-0 #裏
return(men)
}

#－－－－－－ジャンケン(1回発生させる場合)－－－－－－
zyanken<-function()
{
x<-runif(1)
if (x<=1/3) te<-1
else if (x<=2/3) te<-2
else te<-3
return(te)
}

#－－－－－－コイン投げ:coin=1表，0=裏－－－－－－
#－－－－－－n回投げたときに表が何回でるか？－－－－

mysimulation<-function(n)
{
count<-0 #カウンタを初期化する．
for (i in 1:n)
{
x<-runif(1)
if (x<=1/2) coin<-1
else coin<-0
if (coin==1) count<- count+1
}
return(count)

次はモンテカルロシミュレーションの基礎の基礎です．

#－－－－－－コイン投げのサンプル－－－－－－
coin<-function(){
x<-runif(1)
if (x<=1/2) men<-1
else men<-0
return(men)
}

#－－－－－－モンテカルロシミュレーションの単純な例1－－－－－－
#－－－－－－表がでる確率がどうなるか？－－－－－－
mymontecarlo<-function(n){
count<-0
for (i in 1:n){
count<-count+coin()
}
return(count/n)
}

#－－－－－－モンテカルロシミュレーションの単純な例2－－－－－－
#－－－－－－10回コインを投げた平均の平均－－－－－－－－－－
mycointoss<-function(n){
count<-0
for (i in 1:n)
{
x<-coin()
if (x==1) count<-count+1
}
return(count)
}
coin.montecarlo<-function(n){
count<-0
for (i in 1:n)
{
count<-count+mycointoss(10)
}
return(count/n)
}

次は，円周率を求めるサンプルです．こんなことも出来るんですね．
恥ずかしながら，今頃知りました…（汗）．

#－－－－－－モンテカルロシミュレーションの単純な例4－－－－－－
#－－－－－－円周率を求める－－－－－－－－－－－－－
#－－－－－－1/4円を対象としている－－－－－－－－－－－－－
pi.montecarlo<-function(n){
count<-0
for (i in 1 : n){
suna<-runif(2)
if (sqrt(suna[1]^2+suna[2]^2)<1) {
count<-count+1
}
}
return(4*count/n)
}

モンテカルロシミュレーション，サンプルプログラムの続きです．
今度は，以下の二つです．
　・円周率を求める．全部の円を対象として．
　・二つの円（重なり有）の面積の割合を求める．

シミュレーションというのは初めてやったのですが，
新しい使い方を知ることが出来て面白いです♪
（今まで，ぜんぜんやっていないことが，バレマスガ…）

#－－－－－－モンテカルロシミュレーションの単純な例4－－－－－－
#－－－－－－円周率を求める－－－－－－－－－－－－－－－－－－
#－－－－－－全部の円を対象としている－－－－－－－－－－－－－

pi.montecarlo2<-function(n){
count<-0
for (i in 1:n){
suna<-runif(2,min=-1,max=1)
if (sqrt(suna[1]^2+suna[2]^2)<1){
count<-count+1
}
}
return(4*count/n)
}

#－－－－－－モンテカルロシミュレーションの単純な例5－－－－－－
#－－－－－－二つの円（1部重なり有）の面積を求める－－－－－－
circle.montecarlo<-function(n){
count<-0
for (i in 1:n){
suna<-runif(2)
pi=sqrt((suna[1]-0.3)^2 + (suna[2]-0.3)^2)
if (pi<0.3) count<-count+1

else {
pi=sqrt((suna[1]-0.6)^2 + (suna[2]-0.6)^2)
if (pi<0.4) count<-count+1
}
}
return(count/n)
}

次はPlot関数の簡単なSampleを記します．
ＳＡＳよりも，数段，簡単で綺麗なようです．
まあ，そう思ったから覚えようと思ったんだけど…．

#－－－－－－－－Plot関数について－－－－－－－－
plot(x軸のデータ, (y軸のデータ) ,オプション)
plot(関数名，x軸範囲の下限,y軸範囲の下限)

#－－－－－－－－Plot関数，オプション，－－－－－－
# 引数TYPEの説明で，プロットの形式を定義(8種類有)
# 軸の制御の引数： log, xlim, ylim
# タイトル等の設定： main, sub, xlab, ylab, tmag
#－－－－－－－－Plot関数，Sample－－－－－－－－
x<-c(1,2,3,4,5,6,7,8,9,10)
plot(x)

x<-1:10
y<-1:10
plot(x,y)
plot(x,y,xlim=c(10,1))

#－－－－－－－－－－－－－－－－
plot(sin,-pi,2*pi)

#－－－－－－－－－－－－－－－－

gauss.density<-function(x){
1/sqrt(2*pi)*exp(-x^2/2)
}
plot(gauss.density,-3,3)

#－－－－－－－－－－－－－－－－
x<-rnorm(10)
plot(x,type="l")

#－－－－－－－－－－－－－－－－
x<-rnorm(10)
plot(x, ylim=c(-30,30),type="l")

#－－－－－－－－－－－－－－－－
x<-rnorm(10)
plot(x,main="Simple Time Series")

#－－－－－－－－－－－－－－－－
plot(sin,-pi,pi,xlab="x",ylab="y",lty=2)
plot(cos,-pi,pi,add=T)

#－－－－－－－－－－－－－－－－
#－－－－－－Plotを条件分岐する例－－－
x<-runif(100)
y<-runif(100)
plot(x,y,pch=ifelse(y>0.5,1,18))

次は三次元プロットです．今回はあまり使用しないと思うので，さらっと記載します．
（他と同じかも？？）

#－－－－－－－－－－－－－－－－－－－－－－－
#三次元プロット
#persp(x軸データ，y軸データ,z軸データ,col=色,
#theta=横回転角度,phi=縦回転角度,expand=拡大率)
#－－－－－－－－－－－－－－－－－－－－－－－

#二次元正規分布の関数をプロット

x<-seq(-3,3,length=50)
y<-x

gauss3d<-function(x,y){
rho<-0.9
return(1/(2*pi*sqrt(1-rho^2))*exp(-(x^2-2*rho*x*y+y^2)/(2*(1-rho^2))))
}

z<-outer(x,y,theta=30)

persp(x,y,z,theta=30,phi=30,expand=0.5,col="lightblue")

ジャンル別一覧