#===========================
#　VPA　プログラム 
#===========================

VPA<-function(data,M,F.first){
 data<-data
 M<-M
 F.first<-F.first
#まずは普通に計算
	VPA.first<-function(data,M,F.first){
	 C<-data				#漁獲量
	 N<-data				#資源尾数
	 F<-data				#漁獲係数
	 F.first<-F.first			#最新年・最高齢魚のFの設定
	 F[nrow(F),ncol(F)]<-F.first	
	 M<-M					#自然死亡係数

	for (j in 0:(ncol(C)-2)){
		N[nrow(N),ncol(N)-j]<-C[nrow(C),ncol(C)-j]*(F[nrow(F),ncol(F)-j]+M)/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]<-log(N[nrow(N)-i,ncol(N)-i-j]/N[nrow(N)-i+1,ncol(N)-i+1-j])-M
		}
		if ((ncol(C)-j)<=2)break
		F[nrow(F),ncol(F)-j-1]<-F[nrow(F)-1,ncol(F)-j-1]
							#　↑　Fの近似
	}

	for (j in (1:(nrow(C)-1))){
		F[nrow(F)-j,ncol(F)]<-mean(c(F[nrow(F)-j,ncol(F)-1],F[nrow(F)-j,ncol(F)-2],F[nrow(F)-j,ncol(F)-3]))
			# ↑　最新年Fを、3年間の平均値と近似
		N[nrow(N)-j,ncol(N)]<-
			C[nrow(C)-j,ncol(C)]*(F[nrow(F)-j,ncol(F)]+M)/(F[nrow(F)-j,ncol(F)]*(1-exp((-1)*(F[nrow(F)-j,ncol(F)]+M))))
			#　↑　近似した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]<-log(N[nrow(N)-i-j,ncol(N)-i]/N[nrow(N)-i+1-j,ncol(N)-i+1])-M
		}
	}
	A<-list(推定資源尾数=N,推定漁獲係数=F)
	return(A)
 	}
#最適化のための、残差平方和を計算
	VPA2<-function(F.first){
		model<-VPA.first(data,M,F.first)
		(model$推定漁獲係数[nrow(data),ncol(data)]-model$推定漁獲係数[nrow(data)-1,ncol(data)])^2
	}
#準ニュートン法(quasi-Newton method )による最適化
	saiteki<-optim(F.first,VPA2,method="BFGS")

#最適化計算により求められたF.firstを新たに入れ込んで、再計算
	VPA.first(data,M,saiteki$par)
}