幂函数 y = kx^b 当 x<=0 有麻烦，我们可以移 y 轴使 x 都为正数，只考虑 x>0 吧。 被拟合的数据要像拟合曲线，如果不像幂函数，可以考虑其它，例如多项式。 已知点 (xi, yi), i = 1, ..., n，要找一条幂函数曲线 y = kx^b 尽量接近这些点， 其中 k, b 为待定系数。 用最小二乘法找这条曲线： 使这 n 个点的函数差的平方的和 sigma (kxi^b - yi)^2 为极小，就是对这个式子中的 k 和 b 求偏导数 = 0， 得到方程组，解出这个方程组未知数 k, b。但这样得到的方程组是非线性的，解起来麻烦，一 般先变换一下。 对 y = kx^b 两边取对数 ln(y) = ln(k) + bln(x)，令 Y = ln(y)，K = ln(k)，B = b，X = ln(x) '(可见，y, k, x 都不能 <= 0) 得到 Y = K + BX。使 sigma (K + BXi - Yi)^2 为极小， 对这个式子中的 K 和 B 求偏导数，令它们等于 0，得到线性方程组。解出未知数 K, B。 Option Explicit Private Sub Command1_Click() '幂函数 y = kx^b 最小二乘法拟合，编程 Rock www.dayi.net Dim n As Integer, i As Integer Dim x() As Double, y() As Double, x1() As Double, y1() As Double, aa() As Double Dim a(2, 2) As Double, c(2) As Double Dim k As Double, b As Double, yTot As Double, yEve As Double, yErr As Double, r2 As Double n = 5 ReDim x(1 To n), y(1 To n), x1(1 To n), y1(1 To n) '数据点数 x(1) = 1: x(2) = 2: x(3) = 3: x(4) = 4: x(5) = 5 据 x 值 y(1) = 1: y(2) = 1.7: y(3) = 2.1: y(4) = 2.3: y(5) = 2.5 数据 y 值 '例子，数 '例子， a(1, 1) = n For i = 1 To n x1(i) = Log(x(i)) y1(i) = Log(y(i)) a(1, 2) = a(1, 2) + x1(i): c(1) = c(1) + y1(i) 得到的方程组矩阵 '变换 '偏导数为 0 a(2, 2) = a(2, 2) + x1(i) * x1(i): c(2) = c(2) + x1(i) * y1(i) Next 

a(2, 1) = a(1, 2) If Not Solve(a(), c()) Then Print "解方程组出错!" Else k = 2.71828182845905 ^ c(1): b = c(2) 线 y = kx^b For i = 1 To n yEve = yEve + y(i) / n Next For i = 1 To n yTot = yTot + (y(i) - yEve) ^ 2 yErr = yErr + (y(i) - (k * x(i) ^ b)) ^ 2 Next r2 = 1 - yErr / yTot 差 Print Print vbTab & "拟合曲线 y = kx^b = " & r2 End If ReDim aa(1 To n, 1 To 4) For i = 1 To n aa(i, 1) = x(i) aa(i, 2) = y(i) aa(i, 3) = x(i) aa(i, 4) = k * x(i) ^ b Next i MSChart1.ChartData = aa '得到拟合曲 '变换前的拟合误 k = " & k & " b = " & b & " r2 '画图 MSChart1.chartType = VtChChartType2dXY MSChart1.Plot.SeriesCollection(1).ShowLine = False MSChart1.Plot.SeriesCollection(1).SeriesMarker.Show = True MSChart1.Plot.SeriesCollection(1).SeriesMarker.Auto = False MSChart1.Plot.SeriesCollection(1).DataPoints(-1).Marker.Style = VtMarkerStyleCircle MSChart1.Plot.SeriesCollection(1).DataPoints(-1).Marker.Size = 180 End Sub Private Function Solve(a() As Double, b() As Double) As Boolean 

'Gaussian elimination method, coded by www.dayi.net btef (please let this line remain) Dim ii As Integer, i As Integer, j As Integer, n As Integer Dim d1 As Double Dim maxV As Double, maxI As Integer n = UBound(b) For ii = 1 To n - 1 maxV = Abs(a(ii, ii)) maxI = ii For i = ii To n If Abs(a(i, ii)) > maxV Then maxV = Abs(a(i, ii)) maxI = i 'row with max element End If Next i If ii <> maxI Then 'exchange rows d1 = b(maxI) b(maxI) = b(ii) b(ii) = d1 For j = ii To n d1 = a(maxI, j) a(maxI, j) = a(ii, j) a(ii, j) = d1 Next j End If If a(ii, ii) = 0# Then Exit Function End If For i = ii + 1 To n 'elemination d1 = a(i, ii) / a(ii, ii) For j = ii + 1 To n a(i, j) = a(i, j) - a(ii, j) * d1 Next j b(i) = b(i) - b(ii) * d1 Next i Next ii b(n) = b(n) / a(n, n) 'back substitution For i = n - 1 To 1 Step -1 d1 = 0# 

For j = i + 1 To n d1 = d1 + a(i, j) * b(j) Next j b(i) = (b(i) - d1) / a(i, i) Next i Solve = True End Function