Description

Here we test that the extrinsic functions in cppcclib for uni-variate
normal distributions (PDFs and CDFs) work as expected.

Contributor: Steve

Small Model of Type : GAMS

Category : GAMS Test library

Main file : cpplib01.gms

$title Test extrinsic functions in cppcclib (CPPLIB01,SEQ=635)

$onText
Here we test that the extrinsic functions in cppcclib for uni-variate
normal distributions (PDFs and CDFs) work as expected.

Contributor: Steve
$offText
$onDollar

$funcLibIn mvnLib cppcclib

function pdf1 'PDF of univariate normal' / mvnLib.pdfUVN /;
function cdf1 'CDF of univariate normal' / mvnLib.cdfUVN /;

$if not set INFILE $set INFILE uvnFull

alias(*,u1,u2);
set V(u2)
/ x
  f,     f_,     f_a,     f_r
  fx,    fx_,    fx_a,    fx_r
  fxx,   fxx_,   fxx_a,   fxx_r
  g,     g_,     g_a,     g_r
  gx,    gx_,    gx_a,    gx_r
  gxx,   gxx_,   gxx_a,   gxx_r
/;
scalar
 aeps0       'absolute error tolerance: function' / 1e-15 /
 reps0       'relative error tolerance: function' / 1e-15 /
 aeps1       'absolute error tolerance: gradient' / 1e-15 /
 reps1       'relative error tolerance: gradient' / 1e-15 /
 aeps2       'absolute error tolerance: Hessian' / 1e-13 /
 reps2       'relative error tolerance: Hessian' / 1e-13 /
 aepsn       'absolute error tolerance: numeric' / 1e-5 /
 repsn       'relative error tolerance: numeric' / 1e-5 /
 ;
set T(u1);
parameter data(u1,u2);

$gdxIn %INFILE%
$load data
$gdxIn

option T < data;

data(T,'x') = data(T,'x') + eps;
loop {T,
  data(T, 'f'  ) = cdf1.value(data(T,'x'));
  data(T, 'fx' ) = cdf1.grad (data(T,'x'));
  data(T, 'fxx') = cdf1.hess (data(T,'x'));
  data(T, 'g'  ) = pdf1.value(data(T,'x'));
  data(T, 'gx' ) = pdf1.grad (data(T,'x'));
  data(T, 'gx_') = pdf1.gradn(data(T,'x'));
  data(T, 'gxx') = pdf1.hess (data(T,'x'));
  data(T, 'gxx_')= pdf1.hessn(data(T,'x'));
};
data(T,  'fx_'  ) = data(T, 'g');
data(T,  'fxx_' ) = data(T, 'gx');
data(T,  'f_a'  ) = abs(data(T,   'f')-data(T,   'f_'));
data(T,  'fx_a' ) = abs(data(T,  'fx')-data(T,  'fx_'));
data(T,  'fxx_a') = abs(data(T, 'fxx')-data(T, 'fxx_'));
data(T,  'g_a'  ) = abs(data(T,   'g')-data(T,   'g_'));
data(T,  'gx_a' ) = abs(data(T,  'gx')-data(T,  'gx_'));
data(T,  'gxx_a') = abs(data(T, 'gxx')-data(T, 'gxx_'));

set TT(u1);
parameter tmp(u1);

tmp(T) = max[abs(data(T,'f_')),abs(data(T,'f'))];
TT(T) = tmp(T) > 0;
data(T,   'f_r') = INF;
data(T ,  'f_r')$(data(T, 'f_a') eq 0) = 0;
data(TT,  'f_r') = data(TT, 'f_a') / tmp(TT);

tmp(T) = max[abs(data(T,'fx_')),abs(data(T,'fx'))];
TT(T) = tmp(T) > 0;
data(T,  'fx_r') = INF;
data(T , 'fx_r')$(data(T, 'fx_a') eq 0) = 0;
data(TT, 'fx_r') = data(TT, 'fx_a') / tmp(TT);

tmp(T) = max[abs(data(T,'fxx_')),abs(data(T,'fxx'))];
TT(T) = tmp(T) > 0;
data(T,  'fxx_r') = INF;
data(T , 'fxx_r')$(data(T, 'fxx_a') eq 0) = 0;
data(TT, 'fxx_r') = data(TT, 'fxx_a') / tmp(TT);

tmp(T) = max[abs(data(T,'g_')),abs(data(T,'g'))];
TT(T) = tmp(T) > 0;
data(T,   'g_r') = INF;
data(T ,  'g_r')$(data(T,  'g_a') eq 0) = 0;
data(TT,  'g_r') = data(TT, 'g_a') / tmp(TT);

tmp(T) = max[abs(data(T,'gx_')),abs(data(T,'gx'))];
TT(T) = tmp(T) > 0;
data(T,  'gx_r') = INF;
data(T , 'gx_r')$(data(T, 'gx_a') eq 0) = 0;
data(TT, 'gx_r') = data(TT, 'gx_a') / tmp(TT);

tmp(T) = max[abs(data(T,'gxx_')),abs(data(T,'gxx'))];
TT(T) = tmp(T) > 0;
data(T,  'gxx_r') = INF;
data(T , 'gxx_r')$(data(T, 'gxx_a') eq 0) = 0;
data(TT, 'gxx_r') = data(TT, 'gxx_a') / tmp(TT);


set badT(u1);
parameter failures(u1,u2);

badT(T) = ((data(T,   'f_a') > aeps0) and (data(T,   'f_r') > reps0))
      OR  ((data(T,  'fx_a') > aeps1) and (data(T,  'fx_r') > reps1))
      OR  ((data(T, 'fxx_a') > aeps2) and (data(T, 'fxx_r') > reps2))
      OR  ((data(T,   'g_a') > aeps0) and (data(T,   'g_r') > reps0))
      OR  ((data(T,  'gx_a') > aepsn) and (data(T,  'gx_r') > repsn))
      OR  ((data(T, 'gxx_a') > aepsn) and (data(T, 'gxx_r') > repsn))
      ;
failures(badT,V) = data(badT,V);
scalar nTests, nErrors;
nTests = card(T);
nErrors = card(badT);

display 'absolute tolerance, func: ', aeps0;
display 'relative tolerance, func: ', reps0;
display 'absolute tolerance, grad: ', aeps1;
display 'relative tolerance, grad: ', reps1;
display 'absolute tolerance, hess: ', aeps2;
display 'relative tolerance, hess: ', reps2;
display 'failed tests', failures;
* display data;
display 'data points tested: ', nTests;
display '            errors: ', nErrors;
abort$(nErrors) 'There were errors';