Langmuir Model Curve Fit equation

[VB6_to_PBx]

Anyone know how to find the
A, B, C coefficients of this Curve Fit equation ??

Langmuir Model (Extended )
1/(a+b*x^(c-1))

any PureBasic code examples of Curve Fitting ?

my Forum search did not turn up anything useful

[alter Mann]

If I understood that correctly, you want to calculate the coefficients A, B and C from a set of y and x values?
This can be done in the case of non-polynomial functions only with numerical operations.
I've already programmed something similar for another function. If you are interested, I can search the source code. Do you have sample data to try?

[VB6_to_PBx]

If I understood that correctly, you want to calculate the coefficients A, B and C from a set of y and x values?
This can be done in the case of non-polynomial functions only with numerical operations.
I've already programmed something similar for another function. If you are interested, I can search the source code. Do you have sample data to try?

If I understood that correctly, you want to calculate the coefficients A, B and C from a set of y and x values?

Yes , that's correct

X(1) = 60 : Y(1) = 0.999

X(2) = 330 : Y(2) = 2.795

X(3) = 660 : Y(3) = 4.293

X(4) = 1000 : Y(4) = 5.582

X(5) = 1320 : Y(5) = 6.670

using commercial Curve Fitting Software ... i get this for above X and Y Data pairs

Langmuir Model extended : y=1/(a+b*x^(c-1))
Coefficient Data:
a = -0.014900964
b = 11.306459
c = 0.41154387

Code: Select all

Solve for Y
Y = 1 / (a + b * Pow(X,(c - 1)))

However, i need to do this in PureBasic Code instead of using a commercial Curve Fitting Software

any help much appreciated !

[NicTheQuick]

If you change the formula to

Code: Select all

1 / Y = a + b * Pow(X, (c - 1))

and invert all your Y values, it looks like a simple polynomial function which you should be able to fit using the least suqares method.

Edit:
Then there is also this thread about the GNU Scientific Library. Maybe you can find some ideas there.

[alter Mann]

This is my code (also after the suggestion of Nick)

Code: Select all

EnableExplicit

Prototype.d PHI_Fkt (Array c.d(1), x.d)
Prototype.i PHI_Abl (x.d,Array c.d(1), Array a.d(1))

#PB_DBL_EPS        = 2.2204460492503131e-016   ; smallest value such that 1.0+#PB_DBL_EPS != 1.0

Procedure.d PB_NumQuadSum (n.i,Array v.d(1))
Protected dSum.d=0.0;
Protected i.i

For i=0 To n-1 Step 1
  dSum + v(i)*v(i)
Next i
ProcedureReturn dSum

EndProcedure
Procedure.i PbNum_House (iM.i,iN.i,Array dMat.d(2),Array dR.d(1))

  Protected i.i, j.i, k.i
  Protected dFak.d, dRang.d, dAlfa.d, dAkt.d, dEps.d = 2.0 * #PB_DBL_EPS, dTmp.d, dSum.d, dNorm.d, dMaxNorm.d
  Protected Dim dDiag.d(iN)

  If iM < 1 Or iN < 1 Or iM < iN
    ProcedureReturn 1
  EndIf
  
  For i = 0 To iN-1 Step 1 
    dRang = 0.0
    For k = i To iM-1
      dRang + dMat(k,i) * dMat(k,i)
    Next k
    If dRang = 0.0         
      ProcedureReturn 2
    EndIf

    If dMat(i,i) >= 0.0
      dAlfa = Sqr(dRang)
    Else
      dAlfa = - Sqr(dRang)
    EndIf

    dAkt = 1.0 / (dRang + dAlfa * dMat(i,i))
    dMat(i,i) + dAlfa

    dDiag(i) = - dAlfa

    dMaxNorm = 0.0
    For k = i + 1 To iN-1 Step 1
      dFak = 0.0
      dNorm = 0.0
      For j = i To iM-1 Step 1
        dTmp = dMat(j,k)
        dFak + dTmp * dMat(j,i)
        dNorm + dTmp * dTmp
      Next j

      If dNorm > dMaxNorm
        dMaxNorm = dNorm
      EndIf
      dFak * dAkt
      For j = i To iM-1 Step 1
        dMat(j,k) - dFak * dMat(j,i)
      Next j
    Next k

    If Abs(dAlfa) < dEps * Sqr(dMaxNorm)
      ProcedureReturn (3)
    EndIf
    dFak = 0.0
    For j = i To iM-1 Step 1
      dFak + dR(j) * dMat(j,i)
    Next j
    dFak * dAkt
    For j = i To iM-1 Step 1
      dR(j) - dFak * dMat(j,i)
    Next j
  Next i

  For i = iN - 1 To 0 Step -1
    dSum = 0.0
    For k = i + 1 To iN-1 Step 1
      dSum + dMat(i,k) * dR(k)
    Next k
    dR(i) = (dR(i) - dSum) / dDiag(i)
  Next i

  ProcedureReturn 0
EndProcedure

Procedure.i PbNum_NonLinearApprox(m.i,n.i,Array dX.d(1),Array dY.d(1),Array dW.d(1),*PHI.PHI_Fkt,iAblOK.i,*ABL.PHI_Abl,*iMaxIt.integer,dRelEps.d,Array dC.d(1),*dMiQuFe.double)
Protected iIteration.i
Protected i.i, k.i    
Protected Dim dAb.d(m-1,n)
Protected Dim dAbx.d(n)
Protected Dim dB.d(m-1)
Protected Dim dS.d(n)  
Protected Dim dCneu.d(n)
Protected Dim dG.d(m-1) 
Protected dNeuerFehler.d
Protected dDaempfung.d  
Protected dFaktor.d,dCk.d,dHk.d,dZHk.d,dDQ.d

If m <= n Or *iMaxIt\i <= 0 Or dRelEps <= 0.0 : ProcedureReturn 1 : EndIf

For i = 0 To m-1 Step 1                   
  dG(i) = Sqr(dW(i))                     
Next

dFaktor = Pow(#PB_DBL_EPS, 1.0/3.0)

For i=0  To m-1 Step 1
  dB(i) = (dY(i) - *PHI(dC(), dX(i))) * dG(i)
Next i
  
*dMiQuFe\d = PB_NumQuadSum(m, dB())

iIteration = 1
Repeat
  If iAblOK=#True 
    For i=0 To m-1 Step 1                       
      *ABL(dX(i), dC(), dAbx())
      For k=0 To n Step 1
        dAb(i,k) = dAbx(k)
      Next k
    Next i                                               
  Else                         
    For k=0 To n Step 1
      dCk = dC(k)                                
      dHk = dFaktor                               
      If dCk <> 0.0
        dHk * Abs(dCk)
      EndIf
      dZHk = 0.5 / dHk
      For i=0 To m-1 Step 1
        dC(k)     = dCk + dHk
        dDQ       = *PHI(dC(), dX(i))
        dC(k)     = dCk - dHk
        dAb(i,k) = (dDQ - *PHI(dC(), dX(i))) * dZHk
      Next i
      dC(k) = dCk
    Next k
  EndIf
  For i=0 To m-1 Step 1
    For k=0 To n Step 1
      dAb(i,k) * dG(i) 
    Next k
  Next i
  If PbNum_House(m, n + 1, dAb(), dB())
    ProcedureReturn 2 
  EndIf               
  For i=0 To n Step 1 
    dS(i) = dB(i)     
  Next i              
  
  dDaempfung = 1.0
  For k=0 To 10 Step 1 
    For i=0 To n Step 1
      dCneu(i) = dC(i) + dS(i) * dDaempfung   
    Next i
    For i=0 To m-1 Step 1 
      dB(i) = (dY(i) - *PHI(dCneu(), dX(i))) * dG(i)
    Next i
    dNeuerFehler = PB_NumQuadSum(m, dB()) 
    If dNeuerFehler <= *dMiQuFe\d 
      For i=0 To n Step 1       
        dS(i) = dCneu(i) - dC(i)
        dC(i) = dCneu(i)        
      Next i
      Break                     
    EndIf
    dDaempfung * 0.5
  Next k

  If k >= 10                    
    For i=0 To n Step 1
      dC(i) + dS(i)
    Next
    For i=0 To m-1 Step 1       
      dB(i) = (dY(i) - *PHI(dC(), dX(i))) * dG(i)
    Next i
    dNeuerFehler = PB_NumQuadSum(m, dB()) 
  EndIf

  If PB_NumQuadSum(n + 1, dS()) < dRelEps * PB_NumQuadSum(n + 1, dC())
    Break                                  
  EndIf
  If iIteration > *iMaxIt\i
    ProcedureReturn 4
  EndIf

  *dMiQuFe\d = dNeuerFehler
  iIteration + 1
ForEver

*dMiQuFe\d = Sqr(dNeuerFehler)

*iMaxIt\i = iIteration

ProcedureReturn 0
  
EndProcedure

Procedure.d LangmuirFkt (Array dC.d(1), dX.d)
; ***********************************************************
; * the value of a real function with parameters c [i] at x *
; * calculate and return as a function value.               *
; ***********************************************************
  
  ; function changed into 1/y=A+B*x^C
  ProcedureReturn dC(0)+dC(1)*Pow(dX,dC(2))
EndProcedure

Procedure.i LangmuirAbl (dX.d, Array dC.d(1), Array dF.d(1))
; *************************************************************************
; * the partial derivatives f (i) of the above compensation function PHI0 *
; * compute c (i) at the position x according to their parameters.        *
; *************************************************************************
  Protected dPow.d = Pow(dX,dC(2))
  
  dF(0) = 1.0                    ; Derivation of the equation according to the 1st parameter
  dF(1) = dPow                   ; Derivation of the equation according to the 2nd parameter
  dF(2) = dC(1)*dPow*Log(dX)     ; Derivation of the equation according to the 3rd parameter

EndProcedure

Procedure Main()
  Protected iAnz.i,i.i,iIter.i=1000,iFehl.i
  Protected dEps.d=0.0000001,dQEps.d
  Protected Dim dX.d(1)
  Protected Dim dY.d(1)
  Protected Dim dW.d(1)
  Protected Dim dC.d(3)
  
  Read iAnz
  ReDim dX(iAnz)
  ReDim dY(iAnz)
  ReDim dW(iAnz)
  
  For i=0 To iAnz-1 Step 1
    Read.d dX(i)
    Read.d dY(i)
    dY(i) = 1.0/dY(i) ; because of changed function
    Read.d dW(i)
  Next i
  
  dC(0) = 1
  dC(1) = 1
  dC(2) = 1
  
  iFehl = PbNum_NonLinearApprox(iAnz,2,dX(),dY(),dW(),@LangmuirFkt(),#True,@LangmuirAbl(),@iIter,dEps,dC(),@dQeps)
  If iFehl
    MessageRequester("Attention","Error: "+Str(iFehl))
  Else
    MessageRequester("Parameter","A="+StrD(dC(0),5)+"  B="+StrD(dC(1),5)+"  C="+StrD(dC(2)+1.0,5) + "(" + Str(iIter) + ")") 
  EndIf
EndProcedure

; X(1) = 60 : Y(1) = 0.999
; 
; X(2) = 330 : Y(2) = 2.795
; 
; X(3) = 660 : Y(3) = 4.293
; 
; X(4) = 1000 : Y(4) = 5.582
; 
; X(5) = 1320 : Y(5) = 6.670
; 
; using commercial Curve Fitting Software ... i get this For above X And Y Data pairs 
; 
; Langmuir Model extended : y=1/(a+b*x^(c-1))
; Coefficient Data:
; a =	-0.014900964
; b =	11.306459
; c =	0.41154387

Main()

DataSection
  Data.i 5
  Data.d 60.0   ,  0.999 , 1.0
  Data.d 330.0  ,  2.795 , 1.0
  Data.d 660.0  ,  4.293 , 1.0
  Data.d 1000.0 ,  5.582 , 1.0
  Data.d 1320.0 ,  6.670 , 1.0
  
  Data.d -0.014900964 , 11.306459 , 0.41154387
EndDataSection

The result of this code : A=-0.01512 B=11.28654 C=0.41197

[VB6_to_PBx]

alter Mann ,

thanks a bunch for your solution ... it works perfectly !!!

i sent you PM , check it out

Do you need any other Curve Fit Model's equations
i have a ton of Curve Fitting Models + solutions
maybe i might have a particular Curve Fit Model you might be interested in ??

also thank you NicTheQuick for help + links !!!

[alter Mann]

Unfortunately, I deleted all the comments because they were all in German. Therefore here are some explanations:
The code can also be used for other approximations. For this purpose, the functions LangmuirFkt and LangmuirAbl must be replaced and the array dC () must be adapted to the number of unknown parameters. For simple functions any starting values ​​can be used (see dC (0) = 1.0 ... in the source code). With the original function y = 1 / (A + Bx ^ (C-1)) the start values must be close to the solution. By the transformation to Nick's proposal in 1 / y = A + Bx ^ C the calculation was much
more stable and the starting values ​​almost arbitrary.

[NicTheQuick]

I also asked my girlfriend. She said it will work with the least-squares method but you need the non-linear ones, for example the Gauß-Newton algorithm because you can calculate the gradient very nice with it. And in the case it will not converge good enough you can use the gradient free Nelder-Mead method to minimize the quadratic error. That should be very robust.

I have no time to develop a working code but maybe these ideas can help you further.