Levenberg-Marquardt

[james klutho]

Does anyone have any non-linear regression code ported to
PBCC/PBDLL they would be willing to post? Thank you

------------------

[moved from Source Code forum by Administrator]

[John Petty]

James
Here is some code I wrote some years ago.
The code solves multiple linear regressions of up to 8 variables and as a non-linear regression can be converted to a multiple linear regression it also solves multiple or single non -liner regressions.
So if you convert a formula such as:
Y = b0 + b1 x X + b2 x X^2 + b3 x X^3 + b4 x Z + b5 x Z^2
To
Y = b0 + b1 x A0 +b2 x A1 + b3 x A2 + b4 x A3 + b5 x A4
Where
A0 = X, A1 = X^2, A2 = X^3, A3 = Z and A4 = Z^2
Then it can be simply solved as a multiple linear regression of the 5 variables A0 to A4
The following code can handle up to 8 variables of single precision and returns up to 9 b values (b0 to b8) of single precision.

$COMPILE DLL
TYPE X_array
    x(7) AS SINGLE
END TYPE
SUB Dl_Mult_Reg(BYVAL xv AS DWORD, BYVAL yv AS DWORD, np&, nx&, BYVAL bv AS DWORD) EXPORT
'nx&=number of x variables
'np& = number of samples
'bv& points to an array of singles 1 larger than the number of x values
'note the array of x variables must always bi in an array of type X_array
'xvar(np&-1,nx&-1)=x variables
'yvar(np&-1)=y variables
DIM xvar AS X_array pointer
DIM yvar AS SINGLE pointer
DIM b AS SINGLE pointer
REDIM ytot(nx&) AS EXT
REDIM xtot(nx& - 1, nx&) AS EXT    ',0=sum of
REDIM equ(nx&, nx&, nx& + 1) AS EXT
xvar = xv
yvar = yv
b = bv
FOR x& = 0 TO np& - 1
    ytot(0) = ytot(0) + @yvar[x&]
    FOR y& = 1 TO nx&
        ytot(y&) = ytot(y&) + @yvar[x&] * @xvar[x&].x( y& - 1)   'sum of y mult by x's
        xtot(y& - 1, 0) = xtot(y& - 1, 0) + @xvar[x&].x( y& - 1)    'sum of x's
        FOR z& = y& TO nx&
            xtot(y& - 1, z&) = xtot(y& - 1, z&) + @xvar[x&].x( y& - 1) * @xvar[x&].x( z& - 1)
        NEXT
    NEXT
NEXT
equ(0, 0, 0) = ytot(0)
equ(0, 0, 1) = np&
FOR x& = 0 TO nx& - 1
    equ(0, 0, x& + 2) = xtot(x&, 0)
NEXT
FOR y& = 1 TO nx&
    equ(0, y&, 0) = ytot(y&)
    equ(0, y&, 1) = xtot(y& - 1, 0)
    FOR x& = 1 TO nx&
        IF x& < y& THEN
            equ(0, y&, x& + 1) = xtot(x& - 1, y&)
        ELSE
            equ(0, y&, x& + 1) = xtot(y& - 1, x&)
        END IF
    NEXT
NEXT
FOR x& = 0 TO nx& - 1
    l& = nx& + 1 - x&
    FOR y& = 0 TO l& - 2
        FOR z& = 0 TO l& - 1
            equ(x& + 1, y&, z&) = equ(x&, y&, z&) * equ(x&, y& + 1, l&) / 100
            equ(x& + 1, y&, z&) = equ(x& + 1, y&, z&) - equ(x&, y& + 1, z&) * equ(x&, y&, l&) / 100
        NEXT
    NEXT
NEXT
IF equ(nx&, 0, 1) <> 0 THEN @b[0] = equ(nx&, 0, 0) / equ(nx&, 0, 1) ELSE @b[0] = 0
FOR x& = 1 TO nx&
    ltot# = equ(nx& - x&, 0, 0)
    FOR y& = 0 TO x& - 1
        ltot# = ltot# - @b[y&] * equ(nx& - x&, 0, y& + 1)
    NEXT
    IF equ(nx& - x&, 0, x& + 1) <> 0 THEN @b[x&] = ltot# / equ(nx& - x&, 0, x& + 1) ELSE @b[x&] = 0
NEXT
END SUB

The biggest problem is that it is very easy to get overfow or underflow, the program does its best to prevent overflow (the most common) but does not check for either.
typical calling code for say an arrays M, N and O of 100 sets of data for the above equation would be as follows

[CODE]
DIM X(99) as X_array
DIM B(5) as SINGLE
FOR C& = 0 to 99
X(C&).x(0) = N(C&)
X(C&).x(1) = N(c&) * N(C&)
X(C&).x(2) = X(C&).x(1) * N(C&)
X(C&).x(3) = O(C&)
X(C&).x(4) = O(C&) * O(C&)
NEXT
Dl_Mult_Reg X(0), M(0), 100, 5, B(0)

Hope this helps

------------------


[This message has been edited by John Petty (edited August 31, 2001).]

[Peter P Stephensen]

if you are looking for multiple linear regression, have a
look at http://www.powerbasic.com/support/pb...ad.php?t=22848

regards
peter

------------------
[email protected]

[Peter P Stephensen]

If you can use polynomial regression, here is a function

function PolyRegress(x() as double, y() as double, a() as double) as long

    local j as long, p as long, q as long   
    local nObs as long, nPol as long   
 
    nObs = ubound(x())
    if ubound(y()) <> nObs then exit function
 
    nPol = ubound(a())
 
    dim M(nPol,nPol) as local double
    dim b(nPol) as local double
 
    for p = 0 to nPol :	for q = 0 to nPol :	for j = 0 to nObs 
	M(p,q) = M(p,q) + x(j)^(p + q)
	if p = 0 then b(q) = b(q) + y(j) * x(j)^q
    next:next:next
 
    mat M() = inv(M())
    mat a() = M() * b()
 
    function = 1
 
end function

------------------
[email protected]