##########################################################
#####　Richards 関数式への近似によるCuBI height (=hc)の信頼区間を
#####　bootstrapによって求める　                         ver 2010.6.14
##########################################################
#　ReadMe
#     　なお、各行の#の後の部分は注釈とみなされるのでなくてもOKです。
#　準備すること（１）
#　　サンプル用のデータ「TreeDataSample.csv」をダウンロードする。
#　準備すること（２）
#　　Rを起動し、左上の「ファイル」→「ディレクトリの変更...」を選択し、
#　 　ダウンロードした「TreeDataSample.csv」を保存したディレクトリを指定して「OK」選択。
#　　あとは、以下の手順で、指定した部分をRコンソールにコピペするだけ。
#     てっとり早くやりたいなら、とにかく以下全部コピペすればよい

#　手順（１）：データ読み込みなど、下準備
#　　「手順（１）スタート」から「手順（１）エンド」までの行をRコンソールにコピペする

#　手順（２）：bootstrapの計算と計算経過の図示
#　　「手順（２）スタート」から「手順（２）エンド」までの行をRコンソールにコピペする
#　　注１）なお、「手順（２）スタート」のすぐ下の行で、bootstrapを何回やるかを指定する。
#　　　　　デフォルトでは100（100回リサンプリング）になってるが、1000のほうがよい。
#　　注２）再度回数を変えてbootstrapしたいときは、手順（２）からやればよい。

#　手順（３）：bootstrapの計算結果を頻度分布で図示（手順２の結果を使っている）
#　　「手順（３）スタート」から「手順（３）エンド」までの行をRコンソールにコピペする

############　「手順（１）スタート」　################################

# データ読み込み
TestData <-  read.csv("TreeDataSample.csv")

#　読み込んだデータの処理
BA <- (TestData[, "D"]/2/100)^2 *pi     　　　 #胸高直径から個体の胸高断面積 BA (m^2)に換算
TestData2 <- cbind(TestData, BA)　　　　	#　BAを読み込んだデータと結合

#グラフのx軸、y軸の最小値と最大値を生データから決める
xmin <-  0.1*min(TestData2 [,"H"])
xmax <- 1.2 *(max(TestData2 [,"H"]))
ymin <- 0
ymax <- 1.2 * sum(TestData2 [,"BA"])

##　↓Richardsのパラメタ読み込み：bootstrapデータの回帰の初期パラメタとして使う！
##    　　　 生データにRichards関数をあてはめて出した値を初期値に使っている
BB0 <- 0.00930337
AA0 <- 2.167978 
nu0 <-  4.039549 
kk0 <-  2.689324 
hc0 <-  12.10776 

##　サンプリング前にやること：リチャーズ式の非線形回帰を定義関数にした
doit <-  function(cumG, height)
{
  Richards <<- nls(cumG ~   BB + (AA - BB)/(1 + nu * exp(-kk *(height - hc)))^(1/nu) , start = list(BB = BB0, AA = AA0, nu = nu0 , kk = kk0, hc = hc0 ) )
 }
############　「手順（１）エンド」　################################

###########↓ここからブートストラップ###################

############　「手順（２）スタート」　################################

Nboots <- 100  　　　　　 #  bootstrapを何回繰り返すかを指定
                                          # 1000回やると長いので100にしているが、実際は1000にする
bootdata <- c()
namelist <- c("hc") 
# =======
#　　hc3k <- function () {   ###（この行は使ってない）

kkk <- 1

# ↓リサンプリングがNboots（指定したbootstrapの回数）になるまで繰り返す
while(kkk <= Nboots) {  　　　
sampledTree <- array(　, c(length(TestData2[,"H"]),4))
　#　↑個体数分のハコ準備　; [,1]TreeNo    [,2]H    [,3]D   [,4] BA　の順

#  重複を許してデータ数分だけ個体をリサンプリング
TreeNoSample <- sample(length(TestData2[,"TreeNo"]),length(TestData2[,"TreeNo"]), replace = TRUE)
# sample(x, size, ,,, は、ｘから長さsizeでリサンプリングする関数
 # ここでは個体番号（正確にはデータ数分の数値）をサンプリングしている
 
for (nnn in 1: length(TestData2[,1])){
sampledTree[nnn, 1] <- TreeNoSample[nnn]
sampledTree[nnn, 2] <- TestData2[TreeNoSample[nnn], "H"]
sampledTree[nnn, 3] <- TestData2[TreeNoSample[nnn], "D"]
sampledTree[nnn, 4] <- TestData2[TreeNoSample[nnn], "BA"]
}

orderedB <- sampledTree[order(sampledTree[, 2]) ,] 	#樹高の低い順にならびかえ
height <- array(0, dim = c(length(orderedB[,2])) ) 	#樹高を入れるハコを個体数分つくる
cumG <- array(0, dim = c(length(orderedB[,2])) ) 	#BA積算値を入れるハコを個体数分つくる
GGG <- 0

#　個体数分繰り返してｈ－Gデータをつくる。
for (i in 1:length(orderedB[,2])) {
height[i] <- orderedB[i,2]       	#H
GGG <- GGG + orderedB[i,4]  	# orderedB[i,"BA"] は個体のBA
cumG[i] <- GGG　　		# 樹高順にならべた個体のBAを樹高順に積算して並べる（Gの値）
}

#　ｈ－G関係のデータプロット
plot(height, cumG, xlim = c(xmin, xmax), ylim=c(ymin, ymax), mgp = c(2 , 0.6 ,  0 ), xlab = "", ylab = "" )
mtext("Height (m)", side = 1, line = 2)　　　　　　　　　　　#ｘ軸のラベル
mtext(expression(paste("G ( ", m^2, "  )")), side = 2, line = 2) #ｙ軸のラベル
mtext("Bootstrapping for CuBI height (green line) ", side = 3, line = 0.2) # グラフタイトル

#　リチャーズ関数の非線形回帰をやり、エラーが出ないときだけhcを記録
if (class( try (doit(cumG, height), silent = TRUE) ) == "try-error" ) {hc <-  0
} else {try (doit(cumG, height), silent = TRUE) ; hc <-  coef(Richards)["hc"]}
	##   doit(cumG, height)はリチャーズ関数を計算させるルーチン
	###　hc は coef(Richards)["hc"]

# Richards関数にエラーがでなかったときだけ以下を実行
if (hc != 0) {  
curve( coef(Richards)["BB"] + (coef(Richards)["AA"] - coef(Richards)["BB"])/(1 + coef(Richards)["nu"] * exp (-coef(Richards)["kk"] *(x - coef(Richards)["hc"])))^(1/coef(Richards)["nu"]), col="blue", lwd = 2, add = T)
segments(coef(Richards)["hc"], 0, coef(Richards)["hc"], 12, col="green")
text(xmax, 0, paste("N =", kkk), adj = c(1,0)) #　成功したbootstapの回数

bootdata <<- c(bootdata, hc)   	# 個々のhc計算値を記憶
kkk <- kkk +1　　　　　　	#　while 文の回数を数える
}

}　　				#　while(kkk <= Nboots) のおわり　
#　ここですでに繰り返しはおわっている

No.of.data <- length(bootdata)　#　確認のため何回hcを計算させたか調べる
cat(No.of.data)　　　　　　	#　その回数を表示

#　　}　　　			# function hc3k()のおわり（使ってない）

##########################################################
#################　bootstrapping おわり　　  ####################
##########################################################

############　「手順（２）エンド」　################################

#　せっかくシミュレーションしたCuBI Heightの結果（(Hcの値）を保存したければ、ここで
#　write.csv(bootdata, file="RichardsHc.csv")
#   などと入力すると、現在の作業ディレクトリに"RichardsHc.csv"というファイルが作られ
#　1000回分（途中でやめたらそこまでの分）のhcがcsvファイルで保存される。

###########################################################
######　bootstrappingの結果の作図をしたければ、さらに以下を実行する ####
####        1000回シミュレートしたHcの値の頻度分布図が出ます　　########
###########################################################
 
############　「手順（３）スタート」　################################
　
#　上でbootstrapしたhcの値のたまったベクトル「bootdata」を別名で複製
　Hc_values <- bootdata

qqq <- quantile(Hc_values ,p=c(0.025,0.975)) #ｈｃの9５％quantileを取得
LowConfInt <-  qqq[[1]]	#　下から2.5%のhcの値
UppConfInt <-  qqq[[2]]	#　下から97.5%のhcの値を表記用に四捨五入
　LowConfIntPrint <-  sprintf("%1.2f", qqq[[1]])	# 下から 2.5%のhcの値を表記用に四捨五入
 　 UppConfIntPrint <-  sprintf("%1.2f",qqq[[2]])	# 下から97.5%のhcの値を表記用に四捨五入
 　   MeanHc <-  sprintf("%2.2f", mean(Hc_values))	#　hcの平均値を表記用に四捨五入
　   MedianHc <-  sprintf("%2.2f", median(Hc_values))	#　hcの中央値を表記用に四捨五入

#　ヒストグラムのx軸最小値、最大値、頻度階級の間隔（１０等分）の指定
　hxmin <- trunc(min(Hc_values ))
　hxmax <-  ceiling( max(Hc_values ))
　hxinterval <-  (hxmax - hxmin)/10
　histoo <- hist(Hc_values, main = "Bootstrap results for CuBI height", xlab= "", ylab ="", breaks=(seq(hxmin, hxmax ,hxinterval ) ), tcl = -0.5 ,  mgp = c(2 , 0.6 ,  0 ))
 
#　ヒストグラムの軸の説明
mtext("Frequency", side = 2, line = 2)		# 縦軸
mtext("CuBI height  (m) ", side = 1, line = 1.5)             # 横軸

#　その他のデータ表記
mtext( paste("95% quantile = ", LowConfIntPrint, " m  -", UppConfIntPrint  , " m"   ), side = 1, line = 2.5)
mtext( paste("Mean   = ",  MeanHc, " m ;  Median = " ,MedianHc, " m ;  N = ", length(Hc_values) ) , side = 1, line = 3.5)

############　「手順（３）エンド」　################################

###　おしまい　###