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';