DELPHI - How to show chi-square distribution curve

Propability density function (chi-square distribution).

var
  i, iK : integer;
  iX, iY, iGamma, iNumerator, iDenominator : double;
  pSerie : TLineSeries;
  pIntegral : TIntegral;
begin
  Memo1.Clear;

  pSerie := TLineSeries( Chart1.Series[0] );
  pSerie.Clear;

  pIntegral := TIntegral.Create;
  pIntegral.iStep := 0.05;
  pIntegral.SetInterval( 0.001, 20 );

  { -- set param + calc GAMMA function }

  iK := 8;
  iGamma := GetGamma( iK/2 );

  for i := 0 to pIntegral.pX.Count - 1 do
    begin
      iX := pIntegral.pX.GetValue( i );

      try
        iNumerator := power( iX, iK/2-1 ) * exp( ( -1*iX ) /2 );
        iDenominator := power( 2, iK / 2 ) * iGamma;

        if iDenominator <> 0 then
          iY := iNumerator / iDenominator
        else
          iY := 0;
      except
        iY := 0;
      end;
      pIntegral.pY.AddValue( iY );

      { -- add to chart }

      pSerie.AddXY( iX, iY );
    end;
end;

And here is very important Gamma() function with k parameter.

function GammaStirF( X : double ) : double;
var
  y : double;
  w : double;
  v : double;
  stir : double;
begin
  w := 1 / x;
  stir := 7.87311395793093628397E-4;
  stir := -2.29549961613378126380E-4 + w * stir;
  stir := -2.68132617805781232825E-3 + w * stir;
  stir := 3.47222221605458667310E-3 + w * stir;
  stir := 8.33333333333482257126E-2 + w * stir;
  w := 1 + w * stir;
  y := exp(x);
  if x > 143.01608 then
    begin
      v := power( x, 0.5 * x -0.25 );
      y := v *( v / y );
    end
  else
    begin
      y := power( x, x -0.5 ) / y;
    end;
  result := 2.50662827463100050242 * y * w;
end;

function GetGamma( _x : double ) : double;
var
  p : double;
  PP : double;
  q : double;
  QQ : double;
  z : double;
  i : longint;
  SgnGam : double;
begin
  SgnGam := 1;
  q := abs( _x );
  if q > 33.0 then
  begin
    if _x < 0.0 then
    begin
      p := floor(q);
      i := round(p);

      if i mod 2 = 0 then
      begin
        SgnGam := -1;
      end;

      z := q - p;

      if z > 0.5 then
      begin
        p := p+1;
        z := q-p;
      end;

      z := q * sin( pi * z );
      z := Abs(z);
      z := pi / ( z * GammaStirF( q ) );
    end
    else
    begin
      z := GammaStirF( _x );
    end;
    result := SgnGam * z;
    exit;
  end;
  z := 1;
  while _x >= 3 do
  begin
    _x := _x - 1;
    z := z *_x;
  end;
  while _x < 0 do
  begin
    if _x > -0.000000001 then
    begin
      result := z / ( ( 1 + 0.5772156649015329 *_x ) * _x );
      exit;
    end;
    z := z /_x;
    _x := _x + 1;
  end;
  while _x < 2 do
  begin
    if _x < 0.000000001 then
    begin
      result := z / ( ( 1 + 0.5772156649015329 * _x ) *_x );
      exit;
    end;
    z := z /_x;
    _x := _x + 1.0;
  end;
  if _x = 2 then
  begin
    result := z;
    exit;
  end;

  _x := _x - 2.0;
  PP := 1.60119522476751861407E-4;
  PP := 1.19135147006586384913E-3 + _X * PP;
  PP := 1.04213797561761569935E-2 + _X * PP;
  PP := 4.76367800457137231464E-2 + _X * PP;
  PP := 2.07448227648435975150E-1 + _X * PP;
  PP := 4.94214826801497100753E-1 + _X * PP;
  PP := 9.99999999999999996796E-1 + _X * PP;
  QQ := -2.31581873324120129819E-5;
  QQ := 5.39605580493303397842E-4 + _X * QQ;
  QQ := -4.45641913851797240494E-3 + _X * QQ;
  QQ := 1.18139785222060435552E-2 + _X * QQ;
  QQ := 3.58236398605498653373E-2 + _X * QQ;
  QQ := -2.34591795718243348568E-1 + _X * QQ;
  QQ := 7.14304917030273074085E-2 + _X * QQ;
  QQ := 1.00000000000000000320 + _X * QQ;

  result := z * PP / QQ;
  exit;
end;

Output:

k=2

k=8

DELPHI - How to show normal distribution curve (Gaussian)

Typical propability density function (Gaussian distribution).

var
  i : integer;
  iY : double;
  pDN : TDistribution_Normal;
  pSerie : TLineSeries;
begin
  Memo1.Clear;

  /* prepare chart */

  pSerie := TLineSeries( Chart1.Series[0] );
  pSerie.Clear;

  /* prepare distribution data */

  pDN := TDistribution_Normal.Create;
  pDN.Init( 0, 1 );

  for i := 0 to pDN.pIntegral.pY.Count - 1 do
    begin
      { -- get value }

      iY := pDN.pIntegral.pY.GetValue( i );

      pSerie.AddXY( i, iY );
    end;

And here is very important TDistribution_Normal.Init() function with parameters as MEAN and SIGMA, here is generated range MEAN +-3.5 SIGMA.
The sum values under curve is equeal to 1.

function TDistribution_Normal.Init( _iMean, _iSigma : double ) : boolean;
var
  i : integer;
  iFrom, iTo, iValue : double;
begin
  result := false;

  { -- check }

  if _iSigma <= 0 then
    begin
      ShowMessage( 'Distribution_Normal() - standard deviation can`t be 0.' );
      exit;
    end;

  { -- add range -3.5 sigma - Mean - + 3.5 sigma }

  pIntegral.Clear;

  iFrom := _iMean - ( 3.5 * _iSigma );
  iTo := _iMean + ( 3.5 * _iSigma );
  { aproximation - 300 items }
  pIntegral.iStep := ( iTo - iFrom ) / 300;
  pIntegral.SetInterval( iFrom, iTo );

  { -- calc values normal distribution for every X }

  for i := 0 to pIntegral.pX.Count - 1 do
    begin
      iValue := ( 1 / ( _iSigma * sqrt( 2*pi ) ) ) * 
                exp( - ( sqr( pIntegral.pX.GetValue( i ) - _iMean ) ) / 
                ( 2 * sqr( _iSigma ) ) );

      pIntegral.pY.AddValue( iValue );
    end;

  pIntegral.Calc;

  result := true;
end;

Output:

DELPHI - How compute integral - (sin(x))^2 dx from 0 to PI

For integration is used class TIntegral which use for calculation trapezoidal rule.

procedure TForm1.BIntegralClick(Sender: TObject);
var
  i : integer;
  iX, iY : double;
  pIntegral : TIntegral;
  pSerie : TLineSeries;
begin
    pSerie := TLineSeries( Chart1.Series[0] );
    pSerie.Clear;

    { -- prepare integral class - interval 0..PI }

    pIntegral := TIntegral.Create;
    pIntegral.SetInterval( 0, pi );

    for i := 0 to pIntegral.pX.Count - 1 do
      begin
        { calc Y value }
        pIntegral.pY.AddValue( sqr( sin( pIntegral.pX.GetValue( i ) ) ) );

        { -> to graph }
        iX := pIntegral.pX.GetValue( i );
        iY := pIntegral.pY.GetValue( i );
        pSerie.AddXY( iX, iY );
      end;

    ShowMessage( FloatToStr( pIntegral.Calc ) );        // 1.57

    pIntegral.Free;
end;

Here is TIntegral.calc() method:

{ ---------------------------------------------------------------------------
  Method computes integral value (trapezoidal rule).  
  -------------------------------------------------------------------------- }
function TIntegral.Calc( _iCheckSumIntegralValue : double = 999 ) : double;
var
  i, iDeleteXFrom : integer;
  iValue : double;
begin
  result := 0;

  { -- checks }

  if pX.Count < 2 then exit;
  if pX.Count <> pY.Count then exit;

  { -- delete values }

  pIntegral.Clear;
  pCumulation.Clear;

  { -- trapezoidal rule }

  iDeleteXFrom := -1;

  for i := 1 to pX.Count - 1 do
    begin
      { -- compute integral value in interval }

      iValue := (
                  ( pX.GetValueCheckNull( i ) - pX.GetValueCheckNull( i - 1 ) )
                  *
                  ( ( pY.GetValueCheckNull( i ) + pY.GetValueCheckNull( i - 1 ) ) / 2 )
                );

       result := result + iValue;

       { -- check -> is greater than parameter ? }

       if _iCheckSumIntegralValue <> 999 then
       if result > _iCheckSumIntegralValue then
         begin
           { -- subtract }
           result := result - iValue;

           { -- index for removing }
           iDeleteXFrom := i;

           break;
         end;

       { -- record values }

       if bLogIntegral then pIntegral.AddValue( iValue );
       if bLogCumulation then pCumulation.AddValue( result );
     end;

  { -- delete items over test value (if exists) }

  if iDeleteXFrom <> -1 then
    for i := pX.Count - 1 downto iDeleteXFrom do
      begin
        pX.Delete( i );
        pY.Delete( i );
      end;

end;

Output:

And check (same result):

DELPHI - How to make cubic B-Spline interpolation + Cox-De-Boora algorithmus

const giCount = 15;
      giDetail = 10;
...
procedure TForm1.BSpline_Clamped;
var
  j, i, iIndex : integer;
  iNumSegmentsCount : integer;
  iT, i1_T, iB0, iB1, iB2, iB3 : double;
  iNewX, iNewY : double;
  pSerieLine : TLineSeries;
  iValue : extended;
begin

  pSerieLine := TLineSeries( Chart1.Series[1] );
  pSerieLine.Clear;

  { draw count of cubic curves }
  for iIndex := -2 to px.Count - 2 do
    begin

       { every curve generate from I parts }
       for i := 0 to giDetail - 1 do
         begin
           iT := i / (giDetail-1);
           i1_T := 1 - iT;

           iB0 := ( i1_T * i1_T * i1_T ) / 6;
           iB1 := ( ( 3 * iT * iT * iT ) - ( 6 * iT * iT ) + 4 ) / 6;
           iB2 := ( ( -3 * iT * iT * iT ) + ( 3 * iT * iT ) + ( 3 * iT ) + 1 ) / 6;
           iB3 := ( iT * iT * iT ) / 6;

           iNewX := iB0 * GetValueX( iIndex ) +
                    iB1 * GetValueX( iIndex + 1 ) +
                    iB2 * GetValueX( iIndex + 2 ) +
                    iB3 * GetValueX( iIndex + 3 )
                    ;

           iNewY := iB0 * GetValueY( iIndex ) +
                    iB1 * GetValueY( iIndex + 1 ) +
                    iB2 * GetValueY( iIndex + 2 ) +
                    iB3 * GetValueY( iIndex + 3 );

           pSerieLine.AddXY( iNewX, iNewY );
         end;
  end;
  
end;

DELPHI - How to make Chalkin interpolation

procedure TForm1.ChaikinProc( var pNewX, pNewY : TVariable );
var
  iIndex : integer;
  iQx, iQy, iRx, iRy : double;
  iValue1, iValue2 : double;
  pNewListX, pNewListY : TVariable;
begin

  try
    pNewListX := TVariable.Create;
    pNewListY := TVariable.Create;

    { first }
    pNewListX.AddValue( pNewX.GetValue( 0 ) );
    pNewListY.AddValue( pNewY.GetValue( 0 ) );

    for iIndex := 0 to pNewX.Count - 2 do
      begin
        { - calc new Q coord }
        
        iQx := ( 0.75 * pNewX.GetValue( iIndex ) ) + ( 0.25 * pNewX.GetValue( iIndex + 1 ) );
        iQy := ( 0.75 * pNewY.GetValue( iIndex ) ) + ( 0.25 * pNewY.GetValue( iIndex + 1 ) );

        { - calc new R coord }

        iRx := ( 0.25 * pNewX.GetValue( iIndex ) ) + ( 0.75 * pNewX.GetValue( iIndex + 1 ) );
        iRy := ( 0.25 * pNewY.GetValue( iIndex ) ) + ( 0.75 * pNewY.GetValue( iIndex + 1 ) );

        { - add it to values list }

        pNewListX.AddValue( iQx );
        pNewListY.AddValue( iQy );
        pNewListX.AddValue( iRx );
        pNewListY.AddValue( iRy );
      end;

    { last }
    pNewListX.AddValue( pNewX.GetValue( pNewX.Count - 1 ) );
    pNewListY.AddValue( pNewY.GetValue( pNewY.Count - 1 ) );

    { copy }
    pNewListX.CopyTo( pNewX, true );
    pNewListY.CopyTo( pNewY, true );

  finally
    pNewListX.Free;
    pNewListY.Free;
  end;

end;

{ ---------------------------------------------------------------------------
  Chalkin`s curve.
  -------------------------------------------------------------------------- }
procedure TForm1.Chaikin;
var
  i, N : integer;
  pNewX, pNewY : TVariable;
  pSerieLine : TLineSeries;
begin
  pNewX := TVariable.Create;
  pNewY := TVariable.Create;

  for i := 0 to pX.Count - 1 do
    begin
      pX.CopyTo( pNewX );
      pY.CopyTo( pNewY );
    end;

  { get value from scrollbar }
  N := SChalkinFactor.Position;

  { count of smoothing }
  for i := 1 to N do
    ChaikinProc( pNewX, pNewY );

  pSerieLine := TLineSeries( Chart1.Series[1] );
  pSerieLine.Clear;

  for i := 0 to pNewX.Count - 1 do
  begin
    pSerieLine.AddXY( pNewX.GetValue( i ), pNewY.GetValue( i ) );
  end;

  pNewX.Free;
  pNewY.Free;
end;

Smoothing factor = 1

Smoothing factor = 10

DELPHI - How to make Lagrange interpolation

{ ---------------------------------------------------------------------------
  Calc Lagrange polynomial.
  -------------------------------------------------------------------------- }
function TForm1.CalcLagrange( _iX : double ) : double;
var
  i, j : integer;
  iNumerator, iDenominator : double;
begin
  result := 0;

  for i := 1 to pX.Count - 1 do
    begin
      iNumerator := 1;
      iDenominator := 1;

      for j := 1 to pX.Count - 1 do
        if j <> i then
          begin
            iNumerator := iNumerator * ( _iX - pX.GetValue( j ) );
            iDenominator := iDenominator * ( pX.GetValue( i ) - pX.GetValue( j ) );
          end;

      result := result + ( pY.GetValue( i ) * ( iNumerator / iDenominator ) );    
    end;
end;

{ ---------------------------------------------------------------------------
  Lagrange interpolation.
  -------------------------------------------------------------------------- }
procedure TForm1.Lagrange;
var
  j, i, iIndex : integer;
  iNumSegmentsCount : integer;
  iT, i1_T, iB0, iB1, iB2, iB3 : double;
  iNewX, iNewY : double;
  pSerieLine : TLineSeries;
  iValue : extended;
begin

  pSerieLine := TLineSeries( Chart1.Series[1] );
  pSerieLine.Clear;  

  { draw count of cubic curves }
  for iIndex := 0 to px.Count - 2 do
    begin

       { every curve generate from I parts }
       for i := 0 to 10 - 1 do
         begin
           iT := i / (10-1);

           iNewX := pX.GetValue( iIndex ) + 
                    ( (pX.GetValue( iIndex + 1 ) - pX.GetValue( iIndex )) * iT );

           iNewY := CalcLagrange( iNewX );          

           pSerieLine.AddXY( iNewX, iNewY );
         end;

  end;
end;

DELPHI - How make Ferguson cubic interpolation curve

procedure TForm1.Ferguson;
var
  j, i, iIndex : integer;
  iNumSegmentsCount : integer;
  iT, i1_T, iP0x, iP0y, iP1x, iP1y : double;
  iF1, iF2, iF3, iF4 : double;
  iNewX, iNewY : double;
  pSerieLine : TLineSeries;
  iValue : extended;
begin

  pSerieLine := TLineSeries( Chart1.Series[1] );
  pSerieLine.Clear;

  { draw specified count of the cubic curves }
  for iIndex := 0 to px.Count - 2 do
    begin

       { every cube generate from I parts }
       for i := 0 to giDetail - 1 do
         begin
           iT := i / (giDetail-1);

           { -- calculate Hermid polynomials }

           iF1 := ( 2 * power( iT, 3 ) ) - ( 3 * power( iT, 2 ) ) + 1;
           iF2 := ( -2 * power( iT, 3 ) ) + ( 3 * power( iT, 2 ) );
           iF3 := power( iT, 3 ) - ( 2 * power( iT, 2 ) ) + iT;
           iF4 := power( iT, 3 ) - power( iT, 2 );

           { -- and vectors }

           iP0x := abs( pX.GetValue( iIndex + 1 ) - pX.GetValue( iIndex ) );
           iP0y := abs( pY.GetValue( iIndex + 1 ) - pY.GetValue( iIndex ) );

           iP1x := abs( pX.GetValue( iIndex ) - pX.GetValue( iIndex + 1 ) );
           iP1y := abs( pY.GetValue( iIndex ) - pY.GetValue( iIndex + 1 ) );

           iNewX := ( pX.GetValue( iIndex ) * iF1 ) +
                    ( pX.GetValue( iIndex + 1 ) * iF2 ) +
                    ( iP0x * iF3 ) +
                    ( iP1x * iF4 );
           iNewY := ( pY.GetValue( iIndex ) * iF1 ) +
                    ( pY.GetValue( iIndex + 1 ) * iF2 ) +
                    ( iP0y * iF3 ) +
                    ( iP1y * iF4 );

           pSerieLine.AddXY( iNewX, iNewY );
         end;

  end;
end;

DELPHI - How make Catmull-Rom spline (interpolation)

function TForm1.GetCatmull_RomProc( _iT : double; _p0, _p1, _p2, _p3 : double ) : double;
begin
  result := 0.5 * ( ( 2 * _p1 ) +
                    ( ( ( -1 * _p0 ) + _p2 ) * _iT ) +
                    ( ( ( 2 * _p0 ) - ( 5 * _p1 ) + ( 4 * _p2 ) - _p3 ) * _iT * _iT ) +
                    ( ( ( - 1 * _p0 ) + ( 3 * _p1 ) - ( 3 * _p2 ) + _p3 ) * _iT * _iT * _iT )
                  );
end;

{ ---------------------------------------------------------------------------
  Catmull-Rom interpolation curve.
  -------------------------------------------------------------------------- }
procedure TForm1.Catmull_Rom;
var
  j, i, iIndex : integer;
  iNumSegmentsCount : integer;
  iT, i1_T, iP0, iP1, iP2, iP3 : double;
  iNewX, iNewY : double;
  pSerieLine : TLineSeries;
  iValue : extended;
  pNewX, pNewY : TVariable;
begin
  pNewX := TVariable.Create;
  pNewY := TVariable.Create;

  for i := 0 to pX.Count - 1 do
    begin
      pX.CopyTo( pNewX );
      pY.CopyTo( pNewY );
    end;

  { first and last point again }
  pNewX.InsertValue( 0, pNewx.GetValue( 0 ) );
  pNewY.InsertValue( 0, pNewy.GetValue( 0 ) );

  pNewX.InsertValue( pNewX.Count - 1, pNewx.GetValue( pNewX.Count - 1 ) );
  pNewY.InsertValue( pNewY.Count - 1, pNewy.GetValue( pNewY.Count - 1 ) );

  pSerieLine := TLineSeries( Chart1.Series[1] );
  pSerieLine.Clear;

  { draw count of cubic curves }
  for iIndex := 0 to pNewX.Count - 4 do
    begin

       { every curve generate from I parts }
       for i := 0 to giDetail - 1 do
         begin
           iT := i / (giDetail-1);
           i1_T := 1 - iT;

           iNewX := GetCatmull_RomProc( iT, pNewx.GetValue( iIndex ), 
                                        pNewx.GetValue( iIndex + 1 ),
                                        pNewx.GetValue( iIndex + 2 ), 
                                        pNewx.GetValue( iIndex + 3 )
                                      );

           iNewY := GetCatmull_RomProc( iT, pNewy.GetValue( iIndex ), 
                                        pNewy.GetValue( iIndex + 1 ),
                                        pNewy.GetValue( iIndex + 2 ), 
                                        pNewy.GetValue( iIndex + 3 )
                                      );

           pSerieLine.AddXY( iNewX, iNewY );
         end;

  end;

  pNewX.Free;
  pNewY.Free;
end;