# Tuned VPA

Tuned_VPA<-function(data,index,M,F.first){
 data<-data
 index<-index
 M<-M
 F.first<-F.first
#まずは普通に計算
	VPA.first<-function(data,M,F.first){
	 C<-data				#漁獲量
	 N<-data				#資源尾数
	 F<-data				#漁獲係数
	 F[nrow(F)-1,ncol(F)]<-F.first#最高齢一歩手前のFを適当に定める	
	 M<-M					#自然死亡係数
#　最高齢一歩手前のFから、シングルコホート解析する
	N[nrow(N)-1,ncol(N)]<-C[nrow(C)-1,ncol(C)]*( exp(M/2) )/(1-exp(-1*F[nrow(F)-1,ncol(F)]))
							# ↑　まずは最高齢一歩手前魚のNの計算　
	for (i in 1:(nrow(C)-2)){			# ↓　シングルコホート解析
		N[nrow(N)-1-i,ncol(N)-i]<-N[nrow(N)-i,ncol(N)-i+1]*exp(M)+C[nrow(C)-1-i,ncol(C)-i]*exp(M/2)
	}
	for (i in 1:(nrow(C)-2)){			#　↓　Fの計算
		F[nrow(F)-1-i,ncol(F)-i]<- (-1)*log(1-( C[nrow(F)-1-i,ncol(F)-i]*exp(M/2) )/N[nrow(F)-1-i,ncol(F)-i])
	}
#　次は最高齢Fを近似させてから、横にスライドさせていく
	F[nrow(F),ncol(F)]<-F[nrow(F)-1,ncol(F)]
	for (j in 0:(ncol(C)-2)){
		N[nrow(N),ncol(N)-j]<-C[nrow(C),ncol(C)-j]*( exp(M/2) )/(1-exp(-1*F[nrow(F),ncol(F)-j]))
							# ↑　まずは最高齢魚のNの計算　
		for (i in 1:nrow(C)){			# ↓　シングルコホート解析
			if ((ncol(C)-i-j)<=1)break
			N[nrow(N)-i,ncol(N)-i-j]<-N[nrow(N)-i+1,ncol(N)-i+1-j]*exp(M)+C[nrow(C)-i,ncol(C)-i-j]*exp(M/2)
		}
		for (i in 1:nrow(C)){			#　↓　Fの計算
			if ((ncol(C)-i-j)<=1)break
			F[nrow(F)-i,ncol(F)-i-j]<- (-1)*log(1-( C[nrow(F)-i,ncol(F)-i-j]*exp(M/2) )/N[nrow(F)-i,ncol(F)-i-j])
		}
		if ((ncol(C)-j)<=2)break
		F[nrow(F),ncol(F)-j-1]<-F[nrow(F)-1,ncol(F)-j-1]
							#　↑　Fの近似
	}
#　次は近年のFを近似させていく
	for (j in (2:(nrow(C)-1))){
		F[nrow(F)-j,ncol(F)]<-(F[nrow(F)-j,ncol(F)-1]+F[nrow(F)-j,ncol(F)-2]+F[nrow(F)-j,ncol(F)-3])*F[nrow(F)-1,ncol(F)]/(F[nrow(F)-1,ncol(F)-1]+F[nrow(F)-1,ncol(F)-2]+F[nrow(F)-1,ncol(F)-3])
			# ↑　最新年Fを、最新年最高齢一歩手前のFから近似させる
		N[nrow(N)-j,ncol(N)]<-C[nrow(C)-j,ncol(C)]*( exp(M/2) )/(1-exp(-1*F[nrow(F)-j,ncol(F)]))
			#　↑　近似したFから、最新年のNの推定
		if ((nrow(C)-1-j)==0) break
		for (i in 1:(nrow(C)-1-j)){			# ↓　シングルコホート解析
			N[nrow(N)-j-i,ncol(N)-i]<-N[nrow(N)-j-i+1,ncol(N)-i+1]*exp(M)+C[nrow(C)-i-j,ncol(C)-i]*exp(M/2)
		}
		for (i in 1:(nrow(C)-1-j)){			#　↓　Fの計算
			F[nrow(F)-i-j,ncol(F)-i]<-(-1)*log(1-( C[nrow(F)-i-j,ncol(F)-i]*exp(M/2) )/N[nrow(F)-i-j,ncol(F)-i])
		}
	}
	A<-list(推定資源尾数=N,推定漁獲係数=F)
	return(A)
 	}
#最適化のための、残差平方和を計算
	VPA2<-function(F.first){
	Nsum<-apply(VPA.first(data,M,F.first)$推定資源尾数[,-1],2,sum)
	q<-sum(index*Nsum)/sum(Nsum^2)
	sum((index-q*Nsum)^2)
	}
#最適化
	saiteki<-optimise(VPA2,interval=c(0,2))
#最適化計算により求められたF.firstを新たに入れ込んで、再計算
	VPA.first(data,M,saiteki$minimum)
}