################################################
#南風原（2002）第4章
################################################
#Step1:比率の標本分布 p.93-97
################################################

##(a)標本サイズ=5, 母比率=0.6，正答者数=wを作成
N <- 5
PI <- 0.6
w <- seq(0,5,by=1)	#seq(0,1,by=1)は，c(0,1,2,3,4,5)と同様


##(b)4-2式を計算
	#階乗（!）の計算は，gamma(n+1)で計算できる．
	#例えば，5!=5*4*3*2*1の場合は，gamma(5+1)である．
x <- (gamma(N+1)/(gamma(w+1)*gamma(N-w+1)))*PI^w*(1-PI)^(N-w)


##(c)図4-3を描画
	#barplot()で棒グラフを作成できる
	#ヒストグラムとちがい，度数（や比率）が初めからxに含まれている必要がある
	#names.argは，各々のグラフの名前を付ける
barplot(x,ylab="確率",main="正答者wおよび比率pの標本分布",
	names.arg=c(0,1,2,3,4,5))



################################################
#Step2:標本分布の平均 p.97-100
################################################

##(a)正答者数wの標本分布の平均（4.4,4.5式）を確認
sum(x*w)
N*PI


##(b)比率の標本分布の平均（4.6式）を確認
(N*PI)/N


##(c)確率分布の標準偏差（4.7,4.8式）を確認
sqrt(sum((w-N*PI)^2*x))
sqrt(N*PI*(1-PI))


##(e)比率の標本分布の標準偏差（4.9式）を確認
sqrt(N*PI*(1-PI))/N



################################################
#Step3:標準誤差に基づくサンプルサイズの決定 p.101-2
################################################

##(a)標準誤差の関数の定義（4.9式）
##Standard.error(N=標本サイズ,PI=母比率)で標準誤差を算出する
Standard.error <- function(N,PI){
	sqrt(N*PI*(1-PI))/N		#4.8式
	}
Standard.error(5,0.6)	#N=5, PI=0.6の場合のチェック


##(b)図4-4を描画
PI<-seq(0,1,by=0.01)	#PIを0.01刻みで0〜1まで作成する
plot(Standard.error(1,PI)~PI,type="l",
	xlab="母集団比率",ylab="標準誤差",
	main="標本比率pの標準誤差") 		#N=1の場合
lines(Standard.error(3,PI)~PI)			#N=3の場合
lines(Standard.error(10,PI)~PI)			#N=10の場合
lines(Standard.error(100,PI)~PI)		#N=100の場合


##(c)ある標準誤差以下に抑えたい場合の標本サイズを推定
se <- 0.05
PI <- 0.5
(1/se^2)*PI*(1-PI)


##(d)必要な標本サイズを推定する式は便利なので関数を作っておく
req.sample <- function(se,PI){
	(1/se^2)*PI*(1-PI)
	}
req.sample(0.05,0.5)	#標準誤差=0.05，母比率=0.5の場合の標本サイズ



################################################
#Step4:逸脱行動データにおける比率の変動 p.102-3
################################################

##(a)データの読み込み
setwd("C:/data/haebara/chap02")
data1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap1.csv",header=T)	


##(b)変化量が3点以上のデータを抽出
逸脱男 <- subset(data1,性別=="男" & 変化量>=3,select=c(変化量))
逸脱女 <- subset(data1,性別=="女" & 変化量>=3,select=c(変化量))


##(c)比率の推定
	#nrow()は行の数を数えるために使用する
	#length()はベクトルの数を数えるために使用する
	#行列の場合length()は列の数を数えるので注意!!
x <- nrow(逸脱男)/20
y <- nrow(逸脱女)/20


##(d)標準誤差の推定
sqrt(x*(1-x)/20)
sqrt(y*(1-y)/20)


##(e)先に定義した関数を使ってもOK
Standard.error(20,x)
Standard.error(20,y)



################################################
#Step5:正規分布の確率の求め方 p.106-7
################################################

##(a)pnorm()の使用方法は配布資料を参照すること
0.5-pnorm(1.96,lower.tail=F)		#prob(0 < z < 1.96)
pnorm(1.96,lower.tail=F)		#prob(z > 1.96)
0.5-(0.5-pnorm(1.96,lower.tail=F))	#prob(z > 1.96)
(0.5-pnorm(1.96,lower.tail=F))*2	#prob(|z| < 1.96)
1-(0.5-pnorm(1.96,lower.tail=F))*2	#prob(|z| > 1.96)


##(b)pnorm()の使用方法は配布資料を参照すること
qnorm(0.025,lower.tail=F)		#prob(z > zc)=0.025
pnorm(60,mean=50,sd=10,lower.tail=F)	#prob(x >60, mean=50, sd=10)
pnorm((60-50)/10,lower.tail=F)		#prob(x > 1, mean=0,  sd=1)



################################################
#Step6:正規分布を仮定したときの平均の標本分布 p.109-111
################################################

##(a)for(i in 1:n)で，iが1からn番目まで繰り返す
#---------------例---------------------
# y <- seq(3)
# for(i in 1:3){y[i] <- i+5}
# 	ループ1回目：y[1] <- 1+5
# 	ループ2回目：y[2] <- 2+5
# 	ループ3回目：y[3] <- 3+5
# yには 6,7,8が含まれる
#--------------------------------------
x <- rnorm(10000,mean=50,sd=10)	#rnorm(標本サイズ,平均，標準偏差）で
				#大きさが10000の
				#平均50, 標準偏差10の母集団を作成する
hist(x)				

sample(x,10)			#sample(データ,標本サイズ)でデータから無作為に抽出する
sample(x,10)			#無作為だから，1回目と2回目の値は違う

mean(sample(x,10))		#母集団10000から10人サンプリングした場合の標本平均
mean(sample(x,10))		#当然，標本平均も，サンプリングによるバラツキが出る

x.var <- seq(1000)		#大きさが1000のオブジェクトを作る

for(i in 1:1000){		#forループで，1000回サンプリングをして平均を出した
	x.var[i] <- mean(sample(x,10))
	}

x.var				##平均の標本分布（N=10の場合）
hist(x.var)			##平均の標本分布（図4-6の一部）


##(b)平均50，標準偏差10の正規母集団からサンプリングしたとき，
##その標本平均が1以上離れた値をとる確率（N=25の場合）
pnorm(51,mean=50,sd=10/sqrt(25),lower.tail=F)+
pnorm(49,mean=50,sd=10/sqrt(25))


##(c)その標本平均が1以上離れた値をとる確率（N=400の場合）
pnorm(51,mean=50,sd=10/sqrt(400),lower.tail=F)+
pnorm(49,mean=50,sd=10/sqrt(400))


##(d)逸脱行動データにおける平均の変動
setwd("C:/data/haebara/chap02")
data1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap1.csv",header=T)	

逸脱男 <- subset(data1,性別=="男",select=c(変化量))
逸脱女 <- subset(data1,性別=="女",select=c(変化量))
mean(逸脱男)
mean(逸脱女)
mean(逸脱女)-mean(逸脱男)

sd(逸脱男)/sqrt(nrow(逸脱男))
sd(逸脱女)/sqrt(nrow(逸脱女))



################################################
#Step7:2変量正規分布モデル（やらなくても良い） p.112-114
################################################
x <- seq(-3,3,length=40)	#-3〜3まで，長さが40のオブジェクトを作る
y <- seq(-3,3,length=40)	#-3〜3まで，長さが40のオブジェクトを作る


##(b)確率密度関数（4.17式）を作成
binom <- function(x,y,ro){
	sigma.x<-1
	sigma.y<-1
	mu.x   <-0
	mu.y   <-0
	z.x <- (x-mu.x)/sigma.x
	z.y <- (y-mu.y)/sigma.y
	(1/(2*pi*sigma.x*sigma.y*sqrt(1-ro^2)))*
	exp(-((z.x^2-2*ro*z.x*z.y+z.y^2)/(2*(1-ro^2))))
}


##(c)z軸の大きさ（確率密度）を計算する
z <- outer(x,y,ro=0.1,binom)	#outerは(xとyの)外積を計算する．
				#つまり，z軸の大きさを決定する

##(d)persp(x,y,z)で3次元グラフを描く．オプションを適当に変えれば，みる角度が変化する
par(mfrow=c(2,2))		#グラフィックを4分割
persp(x,y,z,theta=150,phi=30,expand=.5,shade=1,col="lightblue",main="2変量正規分布(ρ=0.1)")


z <- outer(x,y,ro=0.6,binom)
persp(x,y,z,theta=150,phi=30,expand=.5,shade=1,col="lightblue",main="2変量正規分布(ρ=0.6)")

z <- outer(x,y,ro=-0.6,binom)
persp(x,y,z,theta=150,phi=30,expand=.5,shade=1,col="lightblue",main="2変量正規分布(ρ=-0.6)")

z <- outer(x,y,ro=0.90,binom)
persp(x,y,z,theta=150,phi=30,expand=.5,shade=1,col="lightblue",main="2変量正規分布(ρ=0.9)")

graphics.off()		#図を閉じる



################################################
#Step8:相関係数の標本分布 p.114-115
################################################

##(a)任意の相関係数を作成する関数を定義 gendat2(標本サイズ,相関係数)
gendat2 <- function(nc, r)	
{
	# 仮のデータ行列を作る。この時点では変数間の相関は近似的に0である。
	z <- matrix(rnorm(2*nc), ncol=2)
	# 主成分分析を行い，主成分得点を求める。この時点で変数間の相関は完全に0である。
	res <- eigen(r2 <- cor(z))
	coeff <-  solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=T))*res$vectors)
	z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff
	# コレスキー分解の結果をもとに，指定された相関係数行列 を持つように主成分得点を変換する。
	z %*% chol(matrix(c(1, r, r, 1), ncol=2))
}


##(b)大きさが10000，相関係数が0.8の母集団を構成する
x <- gendat2(10000,r=0.8)


##(c)母集団からサンプリングをして標本相関係数を算出する
x <- as.data.frame(x)	#xを行列からデータフレームに変換すると
rownames(x)		#行の名前を参照できる

sample(rownames(x),100)	#10000行を100個無作為抽出する
sample(rownames(x),100)	#無作為だから1回目と2回目の値は違う

x[sample(rownames(x),100),]	 #x[行の名前,列の名前]でサンプリングした行を取る
				 #列を指定しない場合は，全ての列が表示される

cor(x[sample(rownames(x),100),]) #標本相関係数を算出する
cor(x[sample(rownames(x),100),]) #当然，標本相関係数も，サンプリングによるバラツキが出る

cor(x[sample(rownames(x),100),])[1,2] #対角行列(1)が邪魔なので，相関係数だけを取る
cor(x[sample(rownames(x),100),])[2,1] #[1,2]でも[2,1]でも良い


##(d)forループで1000回サンプリングをして相関を出す
x.var <- seq(1000)	#大きさが1000のオブジェクトを作る

for(i in 1:1000){
	x.var[i] <- cor(x[sample(rownames(x),100),])[1,2]
	}

x.var			#相関の標本分布（N=100,ρ=0.8の場合）
hist(x.var)		#図4.8の一部


##(e)相関の標本分布の特徴（4.20式）を確認
##標本サイズが大きいほど，母相関係数が大きいほど
##標準誤差が小さくなることが確認できる
ro <- seq(-1,1,by=.01)
plot( (1-ro^2)/sqrt(20) ~ ro, type="l", 
	xlab="母集団相関係数",ylab="標準誤差")	#標本サイズ20
lines( (1-ro^2)/sqrt(50) ~ ro)			#標本サイズ100
lines( (1-ro^2)/sqrt(100) ~ ro)			#標本サイズ100
lines( (1-ro^2)/sqrt(200) ~ ro)			#標本サイズ200
lines( (1-ro^2)/sqrt(300) ~ ro)			#標本サイズ300


##(f)逸脱行動データにおける相関係数の変動
setwd("C:/data/haebara/chap02")
data1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap1.csv",header=T)	

逸脱男 <- subset(data1,性別=="男",select=c(小6,中2))
cor(逸脱男)[1,2]
(1-cor(逸脱男)[1,2]^2)/sqrt(nrow(逸脱男))	#標準誤差



################################################
#Step9:相関係数に関する確率計算 p.115-17
################################################

##(a)フィッシャーのZ変換を定義
fisher.z <- function(r){
		(1/2)*log((1+r)/(1-r))
		}

##(b)図4-9を描画
x <- seq(-1,1,length=50)
plot(fisher.z(x)~x,type="l",
	xlab="相関係数",ylab="Z",main="フィッシャーのZ変換")


##(c)相関係数をZ変換することで標本分布が正規性に近づくことを確認
x <- gendat2(10000,r=0.8)
x <- as.data.frame(x)
x.var <- seq(1000)
for(i in 1:1000){
	x.var[i] <- cor(x[sample(rownames(x),100),])[1,2]
	x.var[i] <- fisher.z(x.var[i])
	}
hist(x.var)		#正規性に近づく


##(d)ρ=.6の2変量正規母集団からN=100のサンプルが取られるとき，r>.7となる確率
se.z <- function(N)(1/sqrt(N-3))	#4.22式を定義
fisher.z(.6)				#r=.6→Z変換
fisher.z(.7)				#r=.7→Z変換
se.z(100)				#N=100の場合の4.23式
pnorm(fisher.z(.7),mean=fisher.z(.6),sd=se.z(100),lower.tail=F)	#解答



################################################
#Step10:回帰係数の標本分布 p.118-19
################################################

##(a)逸脱行動データにおける相関係数の変動
setwd("C:/data/haebara/chap02")
data1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap1.csv",header=T)	
逸脱男 <- subset(data1,性別=="男",select=c(小6,中2))

##(b)残差の標準偏差を算出
reg1 <- summary(lm(中2~小6,data=逸脱男))	#回帰分析の要約をreg1に保存する
reg1		#標準誤差Residual standard errorを確認する
reg1$sigma	#標準誤差は，回帰分析の要約結果のオブジェクト名$sigmaで取り出せる


##(c)標準偏差関数の定義(ssd)と命名
ssd <- function(x){			#BeginFnc
	sqrt(1/length(x)*sum((x-mean(x))^2))
}	
ssd(逸脱男$小6)


#(d)回帰係数の標準誤差
reg1$sigma/(sqrt(nrow(逸脱男))*ssd(逸脱男$小6))
reg1			#実はStd. Errorの欄に表示されている