// regress function
procedure Rat43(const B: TVec; const X: TVecList; const y: TVec);
begin
    // Rat43 := B[0] / Power((1.0 + Exp(B[1]-B[2]*X)),1/B[3]);

    y.Mul(x[0], -B[2]);
    y.Add(B[1]);
    y.Exp();
    y.Add(1);
    y.Power(1/B[3]);
    y.DivideBy(B[0]);
end;

procedure TfrmNonLinTest.FormCreate(Sender: TObject);
var MtxNonLinReg: TMtxMultiNonLinReg;
begin
  MtxNonLinReg := TMtxMultiNonLinReg.Create;
  try
    // Load data - independent variable
    MtxNonLinReg.X.Add(); //add one variable
    MtxNonLinReg.X[0].SetIt(false,[9.000, 14.000, 21.000, 28.000,
                                42.000, 57.000, 63.000, 70.000,
                                79.000]);
    // Load data - dependent variable
    MtxNonLinReg.Y.SetIt(false,[8.930, 10.800, 18.590, 22.330,
                                39.350, 56.110, 61.730, 64.620,
                                67.080]);
    // Initial estimates for regression coefficients
    MtxNonLinReg.b.SetIt(false,[100,10,1,1]);
    // setup optimization parameters
    MtxNonLinReg.Tolerance := 1.0e-6; // 6 digits should do the trick
    MtxNonLinReg.GradTolerance := 1.0e-3 // 3 digits
    MtxNonLinReg.MaxIterations := 400;
    MtxNonLinReg.RegressFunction := Rat43; // regression function
    // Marquardt method
    MtxNonLinReg.OptMethod := optMarquardt;
    MtxNonLinReg.Recalc;
    // MtxNinLinReg.b now stores calculated regression parameter estimates
  finally
    MtxNonLinReg.Free;
  end;
end;

    #include "MtxExpr.hpp"
#include "StatTools.hpp"
#include "Math387.hpp"

void __fastcall Rat43(TVec* const b, TVecList* const x, TVec * const y)
{
    double* B = b->PValues1D(0);

//      return B[0] / Power((1.0 + Exp(b[1]-b[2]*x)),1/b[3]);

    y->Mul((*x)[0], -B[2]);
    y->Add(B[1]);
    y->Exp();
    y->Add(1);
    y->Power(1/B[3]);
    y->DivideBy(B[0]);
}

void __fastcall Example()
{
  TMtxMultiNonLinReg *nlr = new TMtxMultiNonLinReg(NULL);
  try
  {
      const int xLen = 9;
      double xData[xLen] = {9.000, 14.000, 21.000, 28.000, 42.000, 57.000, 63.000, 70.000, 79.000};

      const int yLen = 9;
      double yData[yLen] = {8.930, 10.800, 18.590, 22.330, 39.350, 56.110, 61.730, 64.620, 67.080};

      // Load data - independent variable
      nlr->X->Add(); //add one variable
      (*(nlr->X))[0]->SetIt(false, xData, xLen - 1);

      // Load data - dependent variable
      nlr->Y->SetIt(false, yData, yLen - 1);

      // Initial estimates for regression coefficients
      nlr->B->SetIt(false,OPENARRAY(double,(100,10,1,1)));

      // setup optimization parameters
      nlr->Tolerance = 1.0e-6; // 6 digits should do the trick
      nlr->GradTolerance = 1.0e-5; // 5 digits
      nlr->MaxIteration = 400;
      nlr->RegressFunction = Rat43; // regression function

      // Marquardt method
      nlr->OptMethod = optMarquardt;
      nlr->Recalc();

      // MtxNonLinReg->B now stores calculated regression parameter estimates
  }
  __finally
  {
      nlr->Free();
  }
}