' 發展單位：風禹科技驗證有限公司
' 撰寫人：鄭子璉(Tzu-Lien, Cheng, 璉璉) ，成大水利博肄，微軟最有價值專家
' Web: http://tlcheng.twbbs.org/TLCheng/ E-Mail: qvb3377@ms5.hinet.net
' --------------------------------------------------------------------------------------

Module modStatistic
	Public Enum DataGroup As Integer
		Sample = 0
		[Default] = Sample
		Population = 1
	End Enum

	Public Function GetMaxPercentError(ByVal arrSimulate() As Double, ByVal arrObserved() As Double, Optional ByVal LowerBound As Integer = 0) As Double
		' 最大誤差百分比
		Dim ub, nb, i As Integer
		Dim maxError, nowError As Double

		ub = UBound(arrSimulate)
		nb = ub - LowerBound + 1

		maxError = 0
		For i = LowerBound To ub
			nowError = Math.Abs(arrObserved(i) - arrSimulate(i)) / arrObserved(i)
			If maxError < nowError Then maxError = nowError
		Next

		Return maxError
	End Function

	Public Function GetCoefEfficiency(ByVal arrSimulate() As Double, ByVal arrObserved() As Double, Optional ByVal LowerBound As Integer = 0) As Double
		' 效率係數, Coefficient of Efficiency, CE
		Dim ub, nb, i As Integer
		Dim oMean, summy, summy1, summy2 As Double

		ub = UBound(arrSimulate)
		nb = ub - LowerBound + 1

		oMean = GetAverage(arrObserved, LowerBound)
		For i = LowerBound To ub
			summy1 += (arrObserved(i) - arrSimulate(i)) ^ 2
			summy2 += (arrObserved(i) - oMean) ^ 2
		Next
		summy = 1 - summy1 / summy2

		Return summy
	End Function

	Public Function GetPeekWeightRMSE(ByVal arrSimulate() As Double, ByVal arrObserved() As Double, Optional ByVal LowerBound As Integer = 0) As Double
		' 洪峰均方根物差
		' Peak-weighted root mean square error. This function is identical to
		'the calibration objective function included in computer program HEC-1(USACE, 1998).
		Dim ub, nb, i As Integer
		Dim oMean, summy As Double

		ub = UBound(arrSimulate)
		nb = ub - LowerBound + 1

		oMean = GetAverage(arrObserved)
		For i = LowerBound To ub
			summy += (arrSimulate(i) - arrObserved(i)) ^ 2 * (arrObserved(i) + oMean) / (2 * oMean)
		Next
		summy = Math.Sqrt(summy / nb)
		Return summy
	End Function

	Public Function GetRMSE(ByVal arrSimulate() As Double, ByVal arrObserved() As Double, Optional ByVal LowerBound As Integer = 0) As Double
		' 均方根誤差
		Dim ub, nb, i As Integer
		Dim oMean, summy As Double

		ub = UBound(arrSimulate)
		nb = ub - LowerBound + 1

		For i = LowerBound To ub
			summy += (arrSimulate(i) - arrObserved(i)) ^ 2
		Next
		summy = Math.Sqrt(summy) / nb
		Return summy
	End Function

	Public Function GetLeastSquare(ByVal arrSimulate() As Double, ByVal arrObserved() As Double, Optional ByVal LowerBound As Integer = 0) As Double
		' 最小平方法
		Dim ub, nb, i As Integer
		Dim oMean, summy As Double

		ub = UBound(arrSimulate)

		For i = LowerBound To ub
			summy += (arrSimulate(i) - arrObserved(i)) ^ 2
		Next
		Return summy
	End Function

	Public Function GetCorrelationCoefficient(ByVal arrSimulate() As Double, ByVal arrObserved() As Double, Optional ByVal LowerBound As Integer = 0) As Double
		' 相關係數
		Dim ub, nb, i As Integer
		Dim ax, ay, sx, sy, mxy, v_sum As Double

		ub = UBound(arrSimulate)
		nb = ub - LowerBound + 1

		ax = GetAverage(arrSimulate)
		ay = GetAverage(arrObserved)
		sx = GetStandardDeviation(arrSimulate, , , ax)
		sy = GetStandardDeviation(arrObserved, , , ay)
		mxy = ax * ay

		For i = LowerBound To ub
			v_sum += arrSimulate(i) * arrObserved(i) - mxy
		Next

		Return v_sum / ((nb - 1) * sx * sy)
	End Function

	Public Function LogRegression(ByVal arrX() As Double, ByVal arrY() As Double, Optional ByVal LogBase As Double = 0, Optional ByVal LowerBound As Integer = 0) As Double()
		' y=b*x^a => ln(y) = a*ln(x) + ln(b)
		Dim ub, i, nb As Integer
		Dim rtnArray() As Double

		ub = UBound(arrX)
		Dim logX(ub), logY(ub) As Double

		For i = LowerBound To ub
			nb += 1
			If LogBase > 0 Then
				logX(nb) = Math.Log(arrX(i), LogBase)
				logY(nb) = Math.Log(arrY(i), LogBase)
			Else ' 以自然數為底
				logX(nb) = Math.Log(arrX(i))
				logY(nb) = Math.Log(arrY(i))
			End If
			If Double.IsNegativeInfinity(logX(nb)) OrElse Double.IsNegativeInfinity(logY(nb)) Then nb -= 1
		Next

		ReDim Preserve logX(nb), logY(nb)
		rtnArray = LinearRegression(logX, logY, LowerBound)
		If LogBase > 0 Then ' 計算 b
			rtnArray(1) = LogBase ^ (rtnArray(1))
		Else
			rtnArray(1) = Math.Exp(rtnArray(1))
		End If

		Return rtnArray
	End Function

	Public Function LinearRegression(ByVal arrX() As Double, ByVal arrY() As Double, Optional ByVal LowerBound As Integer = 0) As Double()
		' y=ax+b
		Dim ub, nb, i As Integer
		Dim ax, ay, a_up, a_down, a, b As Double

		ub = UBound(arrX)

		ax = GetAverage(arrX, LowerBound)
		ay = GetAverage(arrY, LowerBound)

		For i = LowerBound To ub
			a_up += (arrX(i) - ax) * (arrY(i) - ay)
			a_down += (arrX(i) - ax) ^ 2
		Next

		a = a_up / a_down
		b = ay - a * ax

		Return New Double() {a, b}
	End Function

	Public Function GetAverage(ByVal arrValue() As Double, Optional ByVal LowerBound As Integer = 0) As Double
		' 平均
		Dim ub, nb, i As Integer
		Dim v_sum As Double

		ub = UBound(arrValue)
		nb = ub - LowerBound + 1

		For i = LowerBound To ub
			v_sum += arrValue(i)
		Next

		Return v_sum / nb
	End Function

	Public Function GetStandardDeviation(ByVal arrValue() As Double, Optional ByVal DataGroupFlag As DataGroup = DataGroup.Default, Optional ByVal LowerBound As Integer = 0, Optional ByVal AverageOfValue As Object = Nothing) As Double
		Dim ub, nb, i As Integer
		Dim av, v_sum As Double

		ub = UBound(arrValue)
		nb = ub - LowerBound + 1

		If IsNothing(AverageOfValue) Then
			av = GetAverage(arrValue)
		Else
			av = AverageOfValue
		End If

		For i = LowerBound To ub
			v_sum = v_sum + (arrValue(i) - av) ^ 2
		Next

		Select Case DataGroupFlag
			Case DataGroup.Default, DataGroup.Sample
				Return Math.Sqrt(v_sum / (nb - 1))
			Case DataGroup.Population
				Return Math.Sqrt(v_sum / nb)
		End Select
	End Function

End Module