이전 글에서 Gauss Jordan 소거법을 위해 연립방정식의 해를 구했다면, 이제 데이터로부터 행렬식, 즉 연립방정식을 작성해야합니다.
행렬방정식은 (x, y) 데이터로부터 만듭니다. 그런데, 행렬의 각 값들은 x의 거듭제곱을 더해서 만듭니다. 만약, 3차 다항식을 만들려면, x값의 6승까지 계산한 후 다시 이들을 더하는 방식으로 만들어야하는데, 불필요하게 반복작업이 됩니다. 데이터 수가 많아지면 중복계산수를 최소화하는 것이 유리합니다.
아래의 함수는 x, y 데이터셋이 배열 변수에 미리 입력되었다고 가정하고, 이 값들로부터 연립방정식을 풀기 위한 행렬을 만들고, 앞서 작성한 연립방정식의 해를 구한 후, 그 결과를 출력해줍니다.
모든 배열은 2차원 배열로 작성하였고, 배열의 index는 모두 1로 시작하도록 하였습니다. 그 이유는 엑셀에서 범위값을 지정하고 값을 일괄적으로 읽어오면, 값은 (1 to n, 1 to m)의 크기를 갖는 variant 배열이기 때문이며, 또한 계산 결과를 출력할 때에도 이러한 형태의 배열일 때 엑셀 시트에 출력하기 편리하기 때문입니다. 대신, 코드에서 항상 2차원 배열로 되어있다보니 조금 군더더기가 있어 보이긴 합니다.
'-------------------------------------------
Function PolynomialFit(iX(), iY(), iOrder As Long) As Variant()
Dim i As Long, j As Long, k As Long, tMatX(), tMat(), tX(), tY(), tA()
Dim tVal As Double, tAve As Double, n As Long, tN As Long
'입력된 데이터의 x, y 데이터 갯수를 변수에 담아둡니다. 데이터 갯수를 확인하기 위해 매번Ubound 함수를 계속 쓰는 것보다 코드가 단순해집니다. 또한, 취급하는 행렬의 크기가 (fitting 차수+1)가 되므로 이 값도 변수에 입력해둡니다.
n = UBound(iX, 1)
tN = iOrder + 1
'x값의 평균 대비, 변동폭이 매우 작은 경우를 대비하여 x값의 평균을 구한 후, 평균을 0점이 되도록 평행이동시킵니다.
tAve = 0
For i = 1 To n
tAve = tAve + iX(i, 1)
Next
tAve = tAve / n
ReDim tX(1 To n, 1 To 1)
For i = 1 To n
tX(i, 1) = iX(i, 1) - tAve
Next
'tMatX는 x_i^0승, x_i^1승, x_i^2승... fitting 차수까지의 각 x에 대한 거듭제곱값을 입력해둡니다.
ReDim tMatX(1 To tN, 1 To n)
For j = 1 To n
tMatX(1, j) = 1
For i = 2 To tN
tMatX(i, j) = tMatX(i - 1, j) * tX(j, 1)
Next
Next
'tMat은 tMatX 배열을 2번 곱해서, tN x tN 행렬을 만듭니다. 각 행마다 거듭제곱을 한 후 더하는 것보다 거듭제곱값을 미리 다 계산해두었으니, 중복 계산은 조금 줄어들 수 있습니다.
ReDim tMat(1 To tN, 1 To tN)
For i = 1 To tN
For j = 1 To tN
tVal = 0
For k = 1 To n
tVal = tVal + tMatX(i, k) * tMatX(j, k)
Next
tMat(i, j) = tVal
Next
Next
'x값의 거듭제곱과 y값의 곱 (행렬방정식의 우측 벡터)을 동일한 방식으로 구해줍니다.
ReDim tY(1 To tN, 1 To 1)
For i = 1 To tN
tVal = 0
For j = 1 To n
tVal = tVal + tMatX(i, j) * iY(j, 1)
Next
tY(i, 1) = tVal
Next
'tMat * tA = tY 라는 행렬방정식을 풀기 위해 Gauss-Jordan 소거법을 이용합니다. 정수로 된 특이한 경우가 아니라면, 실수값을 갖는 여러개의 실측 데이터에서는 tMat의 행렬식이 0이 되지 않기 때문에 반드시 해가 존재합니다. 그러나 혹여라도 그렇지 않은 경우도 있으니 해가 있는지 확인합니다.
If Solve_GaussJordan(tMat, tY, tA) Then
'만약 해가 존재한다면, 행렬방정식을 풀면, tA, 즉 다항식으로 fitting했을 때, 계수 a0, a1, a2,... an을 구할 수 있습니다. 그러나, 이 결과는 모든 x가 tAve만큼 평행이동했을 때 fitting된 결과, 즉 y = a0+a1.(x-tAve)+a2.(x-tAve)^2+...+an.(x-tAve)^n 으로 fitting 된 것이므로, 이 다항식을 모두 전개해주어야 합니다.
'다항식 전개는 파스칼 수열, 또는 파스칼의 삼각형으로부터 구할 수 있습니다.
ReDim tMat(1 To tN, 1 To tN)
For i = 1 To tN
tMat(1, i) = 1
tMat(i, i) = 1
Next
For i = 2 To tN
For j = i + 1 To tN
tMat(i, j) = tMat(i - 1, j - 1) + tMat(i, j - 1)
Next
Next
'파스칼 수열을 구했다면, 모든 다항식을 전개했을 때, x^0, x^1, x^2, ... , x^n의 계수를 구할 수 있습니다. tA 배열을 사용해야하기 때문에, 더이상 사용하지 않는 tY를 재사용하여 다항식이 전개했을 때의 계수들을 저장하여 반환합니다.
tVal = 1
For j = 2 To tN
tVal = tVal * (-tAve)
For i = 1 To tN - j + 1
tMat(i, i + j - 1) = tMat(i, i + j - 1) * tVal
Next
Next
ReDim tY(1 To tN, 1 To 1)
For i = 1 To tN
For j = 1 To tN
tY(i, 1) = tY(i, 1) + tMat(i, j) * tA(j, 1)
Next
Next
Else
'만약 오류가 있다면, 빈 배열을 반환합니다.
ReDim tY(1 To tN, 1 To 1)
End If
PolynomialFit = tY
Erase tMatX, tMat, tX, tY, tA
End Function
'-------------------------------------------
Sub Test()
Dim tX(), tY(), tA()
tX = Selection.Columns(1).Value
tY = Selection.Columns(2).Value
tA = PolynomialFit(tX, tY, 3)
WriteMat tA, ActiveSheet, Cells(1, 4)
End Sub
'-------------------------------------------
엑셀의 A열에는 x값을, B열에는 3차식에 해당하는 데이터가 입력되어 있습니다. 차트를 그려서 3차 다항식으로 추세선을 넣게 되면, 아래 그림과 같이 차트에 표시가 됩니다.
위에서 작성한 함수가 제대로 작동하는지 확인하기 위해, Test 매크로를 실행합니다. A열을 x로, B열을 y로 입력되며, fitting된 결과를 D열에 출력하도록 합니다. 출력된 결과는 a0, a1, a2, a3값으로 차트의 추세선식과 동일한 것을 알 수 있습니다. 위의 예제에서 WriteMat은 배열변수를 시트에 쓰는 함수이며, 이전 글에 작성된 매크로입니다. (바로가기 : https://jhjungx.blogspot.com/2021/04/3.html)
여러개의 열을 선택해서 x, y 열로 지정하고, 위의 함수를 반복해서 적용시키게 되면, 일일이 차트를 그리고, 추세선을 추가할 필요없이 한꺼번에 fitting이 가능해집니다. x, y값의 형식이 2차원 배열 형식으로 해두는 이유가 엑셀에서 Range.Value를 써서 한꺼번에 읽어온 데이터를 별도의 형식으로 변환하지 않고 바로 사용할 수 있도록 하기 위함입니다.
다음 글에서는 x, y열을 반복적으로 지정하는 매크로를 작성하고, fitting한 결과를 출력하는 예를 들어보도록 하겠습니다.
댓글 없음:
댓글 쓰기
의견이나 질문이 있으신 분은 언제든지 댓글을 달아주세요~