!   Fit a discrete distribution to             (FitDiscDist.lng)
  a given target/theoretical distribution.
  Match the theoretical pdf with the discrete
  as closely as possible, subject to exactly matching
  the mean and SD;
! Keywords: Fit distribution, Distribution fit, Probability computation,
            Stochastic;
SETS:
 point : p, pt, errup, errdn;
 distset / N, G, T/; ! Names for target distributions supported;
 distn( distset);
ENDSETS
DATA:
  mu = 10;   ! Target mean of discrete distn;
  sd = 3;    ! Target SD (for Normal & Gamma);
  lo = 1;    ! Lowest value in discrete distn;
  hi = 20;   ! Highest value in discrete distn;
  distn = N; ! Choose a target from N(ormal), G(amma), or T(riangle);
  point = 1..30; ! Workspace;
ENDDATA

SUBMODEL fitdist:
!   Choose the fitted probabilities p( i) to closely 
  match the given target probabilities pt( i), and
  exactly match the target mu and target sd;
!   Minimize the total error of fit;
 min = errtot;
!   Match the mean;
 mu = @SUM( point( i) : i* p( i));
!   Match the sd;
 sd*sd = @SUM( point( i): p(i)*( i - mu)^2);

! Compute the relative errors in matching the theoretical pdf;
 @FOR( point( i):
   ( errup( i) - errdn( i))* pt( i) = p( i) - pt( i);
     );
! Compute total sum of absolute relative errors;
  errtot = @SUM( point( i): errup( i) + errdn( i));
ENDSUBMODEL


procedure donormal:
! Generate a discretized version of the Normal, that
 may however not quite match the mean and SD;
! Take care of low bin;
 prev =  @PNORMCDF( mu, sd, lo + 0.5);
 pt( lo) = prev;
! The interior bins;
 @FOR( point( i) | lo #lt# i #and# i #lt# hi:
   ptnu = @PNORMCDF( mu, sd, i + 0.5); 
   pt( i) = ptnu - prev;
   prev = ptnu;
     );
! Take care of high bin;
  pt( hi) =  1 - prev;
ENDprocedure

procedure dogamma:
! Generate a discretized version of the Gamma, that
 may however not quite match the mean and SD;
! We need shape*scale = mu;
!         shape*scale^2 = sd^2;
! Solving ;
  scale = sd* sd/ mu;
  shape = mu/ scale; 
! Take care of low bin;
  prev = @PGAMMCDF( scale, shape, lo + 0.5);
  pt( lo) = prev;
! The interior bins;
  @FOR( point( i) | lo #lt# i #and# i #lt# hi:
    ptnu = @PGAMMCDF( scale, shape, i + 0.5) ;
    pt( i) = ptnu - prev;
    prev = ptnu;
     );
! Take care of high bin;
  pt( hi) =  1 - prev;
ENDprocedure

procedure dotriangular:
! Generate a triangular distribution in [lo, hi];
! We need: lo + hi + mode = 3*mu,
    lo^2 + mode^2 + hi^2 - lo*mode - lo*hi - mode* hi = 18*sd*sd; 
! Solving ;
  mode = 3*mu - lo - hi;
  sd = ((lo^2 + mode^2 + hi^2 - lo* mode - lo* hi - mode* hi)/18)^0.5; 
! Take care of low bin;
  prev = @PTRIACDF( lo, hi, mode, lo + 0.5);
  pt( lo) = prev;
! The interior bins;
  @FOR( point( i) | lo #lt# i #and# i #lt# hi:
    ptnu = @PTRIACDF( lo, hi, mode, i + 0.5);
    pt( i) = ptnu - prev;
    prev = ptnu;
     );
! Take care of high bin;
  pt( hi) =  1 - prev;
ENDprocedure

CALC:
! Generate the appropriate target distribution;
 @FOR( distn( d) | d #eq# 1: donormal);
 @FOR( distn( d) | d #eq# 2: dogamma);
 @FOR( distn( d) | d #eq# 3: dotriangular);

! Take low and high points out of the problem;
 @FOR( point( i) | i #lt# lo #or# i #gt# hi:
   p( i) = 0;
   pt( i) = 0;
     );
! @gen( fitdist);  @SOLVE( fitdist);
 @CHARTBAR( 'Fitting a discrete distribution', 'x-axis', 'Probability', 'Theoretical', pt, 'Fitted', p);
ENDCALC