##########################################################
###　Richards 関数式をｈーG曲線にフィットさせる　　　　　　　　　　　  ##
# 　　　　　　　　　　　　　　　　　　　　　　                ver 2010.6.11
##########################################################
#　準備すること（１）
#　サンプル用のデータ「TreeDataSample.csv」をダウンロードする。
#　　データの調整を自分でやりたいときだけ、次のようにしてｃｓｖファイルを準備する
#　　・個体番号、樹高（ｍ）、胸高直径（ｃｍ）のデータを入力し、CSVファイルで保存する。
#　　・データを「TreeDataSample.csv」というファイル名で保存する
#　　・このファイルの1行目の１～3列目に以下の変数名をいれる
#  　 ・TreeNo    H        D  
#　　・２行目以下にはデータが入れてある。欠損データは除いておく
#　　・個体番号は番号順になっている必要はないが、同じ番号はないようにする。
#　準備すること（２）
#　　Rを起動し、左上の「ファイル」→「ディレクトリの変更...」を選択し、
#　 　ダウンロードした「TreeDataSample.csv」を保存したディレクトリを指定して「OK」選択。
#　　あとは、以下の部分をすべてRコンソールにコピペするだけ。
#     　なお、各行の#の後の部分は注釈とみなされるのでなくてもOKです。

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

#　読み込んだデータの処理
BA <- (TestData[, "D"]/2/100)^2 *pi     　　　#胸高直径から個体の胸高断面積 BA (m^2)に換算
TestData2 <- cbind(TestData, BA)　　　　　　	        #　BAを読み込んだデータと結合
orderedB <-TestData2[order(TestData[, "H"]) ,]          # 樹高の低い順にデータをならびかえ
Height <- array(0, dim = c(length(orderedB[,"H"])) )  # 樹高を入れるハコを個体数分つくる
cumG <- array(0, dim = c(length(orderedB[,"BA"])) ) # BA積算値(=G)を入れるハコを個体数分つくる
GGG <- 0

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

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

#　ｈ－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) #ｙ軸のラベル

##　リチャーズ関数の初期値の設定
#　       自分で準備したデータを使うときは、以下のnu0, kk0, hc0を指定しなおしてください！
BB0 <-  min(orderedB[,"BA"])　　　	#　BB0はBAの最小値に設定
AA0 <-  sum(orderedB[,"BA"])*1.1　	#　AA0はBAの総和(=Gの最大値）の1.1倍に設定
nu0 <-  3 　　　　　　　　　　　　	#3.3　この設定は難しい
kk0 <-  2　　　　　　　　　　　　 	#2.1　この設定も難しい
hc0 <-  12 　　　　　　　　　　　	#ｈ－G関係をみて、ほぼ変曲点の高さぐらいにhc0を設定

##　リチャーズ関数の非線形回帰
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 ) )

##　回帰パラメタを使ってリチャーズ曲線を描く
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", add = T, lwd = 2)

##　回帰パラメタを使ってCuBI Heightのところに線を入れる曲線を描く
segments(coef(Richards)["hc"], 0, coef(Richards)["hc"], ymax, col="green")

# 計算したリチャーズパラメタを表示用に四捨五入
HcValue <- sprintf("%3.2f", coef(Richards)["hc"]) # CuBI Heightを小数点以下2桁にする
Avalue <- sprintf("%3.3f", coef(Richards)["AA"])
Bvalue <- sprintf("%2.5f", coef(Richards)["BB"])
Nuvalue <- sprintf("%3.3f", coef(Richards)["nu"])
Kvalue <- sprintf("%3.3f", coef(Richards)["kk"])

#CuBI Heightを図に記入
text(HcValue, 0, " CuBI Height ", adj =c(0,0))
text(HcValue, 0, paste(" =  ", HcValue, " m ") , adj =c(0,1))

# 計算したリチャーズパラメタを図中に表示
text(xmax * 0.1, ymax * 0.9, paste(" A = ", Avalue), adj =c(0,0))
text(xmax * 0.1, ymax * 0.8, paste(" B = ", Bvalue), adj =c(0,0))
text(xmax * 0.1, ymax * 0.7, paste(" k = ", Kvalue), adj =c(0,0))
text(xmax * 0.1, ymax * 0.6, paste(" nu = ", Nuvalue), adj =c(0,0))
text(xmax * 0.1, ymax * 0.5, paste(" Hc =", HcValue), adj =c(0,0))

#　おわり

#　AICや対数尤度をコンソールに表示したければ以下のように入力する
AIC(Richards)
logLik(Richards)[1]

#　CuBI Heightをコンソールに表示したければ以下のように入力する
coef(Richards)["hc"]