Announcement

Collapse

Forum Guidelines

This forum is for finished source code that is working properly. If you have questions about this or any other source code, please post it in one of the Discussion Forums, not here.
See more
See less

Logistic / Sigmoid / S Curve fitting

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Logistic / Sigmoid / S Curve fitting

    Based on the discussion in this thread:
    https://forum.powerbasic.com/forum/u...ata#post810968
    (Please make any comments there)

    This application fits a logistic curve to a set of data points,
    Click image for larger version  Name:	LogisticCurve.jpg Views:	0 Size:	24.0 KB ID:	811065

    It displays the best fit function and then displays how good the fit is graphically:

    Click image for larger version  Name:	BestFitMsgbox.jpg Views:	0 Size:	19.5 KB ID:	811063

    Click image for larger version  Name:	BestFitGraphic.jpg Views:	0 Size:	44.5 KB ID:	811064


    Attached zip file contains sample tab delimited data "Testdata.tab"

    '
    Code:
    #COMPILE EXE
    #DIM ALL
    FUNCTION PBMAIN() AS LONG
        LOCAL x AS LONG 'loopcount
        LOCAL lngDatapoint AS LONG
        LOCAL strFile, strT AS STRING
        LOCAL MaxY AS EXT ' Max Y value in data ' to check whether inflection point in  data
        LOCAL L,X0,K AS EXT ''formula parameters
        LOCAL DiffSqr,DiffSqrTot,rSqr,rStore,KStore AS EXT ' formula paramters
        LOCAL y1,y2,x1,x2   AS EXT 'used to find inflection point in data
    
        DIM arrX(1 TO 100) AS EXT ' Use separate arrays for ease of sort
        DIM arrY(1 TO 100) AS EXT ' If initial data is not sorted
    
        DISPLAY OPENFILE , , , "Select data file", EXE.PATH$, CHR$("All Files", 0, "*.*", 0),"Testdata.tab","" , 0 TO strFile
        IF NOT ISFILE(strFIle) THEN
            ? "File Not Found - Exiting"
            EXIT FUNCTION
        END IF
        OPEN strFile FOR INPUT AS #1
        WHILE NOT EOF(#1)
            INCR lngDatapoint
            IF UBOUND(arrX()) < lngDatapoint THEN REDIM PRESERVE arrX(lngDataPoint + 100):REDIM PRESERVE arrX(lngDataPoint + 100)
            LINE INPUT #1, strT
            arrx(lngDatapoint) = VAL(PARSE$(strT,$TAB,1))
            arrY(lngDatapoint) = VAL(PARSE$(strT,$TAB,2))
            IF arrY(lngDataPoint) > MaxY THEN MaxY = arrY(lngDatapoint)
        WEND
        CLOSE #1
    
        REDIM PRESERVE arrX(1 TO lngDatapoint)
        REDIM PRESERVE arrY(1 TO lngDatapoint)
        ARRAY SORT ArrX(), TAGARRAY ArrY()
        'Get limits and inflexion point
        L = VAL(INPUTBOX$("Enter upper limit fo S curve",,STR$(MaxY)))
    
        IF MaxY < L/2 THEN ' no inflection point found
            X0 = VAL(INPUTBOX$("Enter Inflection point (Value of X where Y = " & STR$(L/2) & ")",,""))
        ELSE
            FOR x = 1 TO UBOUND(arrX())
                y2 = arrY(x): x2 = arrX(x)
                IF y2 >= L/2 THEN EXIT FOR
                y1 = y2 :x1= x2
            NEXT
            X0 = x1+ (x2 - x1) * (L/2 - y1)/(y2-y1) 'inflection point.
        END IF
         rSTore = L^2 ' Max R Squared.
    
        FOR k = 0.0001 TO 0.9999 STEP 0.0001
           DiffSqrTot = 0
           FOR x = 1 TO UBOUND(arrX())
               DiffSqr  = (arrY(x) - L/(1+EXP(-(K*(arrX(x)-X0)))))^2
               DiffSqrTot += DiffSqr
           NEXT
           rSqr = DiffSqrTot/UBOUND(arrX())
           IF rSqr < rStore THEN rSTore = rSqr: KSTore = k
       NEXT
       ? STR$(UBOUND(arrX())) & " data points read from " & $LF & strFile & $LF & $LF & _
             "Best Fit: " & $LF & _
             "Y = " & FORMAT$(L) & "/(1 + e^(-" & FORMAT$(kstore,"0.00000") & " * (x - " & FORMAT$(X0,"0.00000") & "))" & $LF & $LF & _
             "R^2 = " & FORMAT$(rStore,"0.00000"),,"Logistic Curve Fit"
    
          Plotit ArrX(),ArrY(),L,kstore,X0,strFile
    END FUNCTION
    
    FUNCTION Plotit(ArrX() AS EXT, arrY() AS EXT, l AS EXT,k AS EXT,x0 AS EXT, strF AS STRING) AS LONG
        LOCAL hWin AS DWORD
        LOCAL strCaption AS STRING
        LOCAL x AS DOUBLE
        strCaption = strF & " (Blue) v Y = " & FORMAT$(L) & "/(1 + e^(-" & FORMAT$(k,"0.00000") & " * (x - " & FORMAT$(X0,"0.00000") & ")) (Red)"
        GRAPHIC WINDOW NEW  strCaption, 50, 50, 800, 300  TO hWin
        GRAPHIC SCALE (0,L) - (ArrX(UBOUND(arrx())),0)
        GRAPHIC COLOR %BLUE,-1
        GRAPHIC WIDTH 3
        GRAPHIC SET POS (arrx(1),arry(1))
        FOR x = 2 TO UBOUND(arrx())
            GRAPHIC LINE - (arrX(x),arrY(x))
        NEXT
        GRAPHIC COLOR %RED,-1
        GRAPHIC WIDTH 3
        GRAPHIC SET POS (arrX(1),L/(1+EXP(-(K*(arrX(1)-X0)))))
        FOR x = 1 TO UBOUND(arrx())
            GRAPHIC LINE - (ArrX(x),L/(1+EXP(-(K*(arrX(x)-X0)))))
        NEXT
         GRAPHIC WAITKEY$
        GRAPHIC DETACH
        GRAPHIC WINDOW END hWin
    END FUNCTION
    '

    Attached Files

  • #2
    Stuart

    Forgive my ignorance

    But if I do this, how do I determine the n+1 figure?

    Thanks

    Kerry
    [I]I made a coding error once - but fortunately I fixed it before anyone noticed[/I]
    Kerry Farmer

    Comment


    • #3
      Originally posted by Kerry Farmer View Post
      Stuart

      Forgive my ignorance

      But if I do this, how do I determine the n+1 figure?

      Thanks

      Kerry
      Taken to the original thread so as not to contaminate the source code forum

      Comment

      Working...
      X