Description
This test makes sure that the random number generators, the probability density functions, the cumulative distribution functions and the inverse cumulative distribution functions from stodclib work as expected. Contributor: L. Westermann
Small Model of Type : GAMS
Category : GAMS Test library
Main file : stolib01.gms
$title Test extrinsic functions in stodclib (STOLIB01,SEQ=520)
$onText
This test makes sure that the random number generators, the probability density
functions, the cumulative distribution functions and the inverse cumulative
distribution functions from stodclib work as expected.
Contributor: L. Westermann
$offText
$funcLibIn stolib stodclib
function dseed /stolib.SetSeed /;
function dunif /stolib.Duniform /;
function pdfunif /stolib.pdfuniform /;
function cdfunif /stolib.cdfuniform /;
function icdfunif /stolib.icdfuniform /;
function dnorm /stolib.Dnormal /;
function pdfnorm /stolib.pdfnormal /;
function cdfnorm /stolib.cdfnormal /;
function icdfnorm /stolib.icdfnormal /;
function dinvgaus /stolib.Dinvgaussian /;
function pdfinvgaus /stolib.pdfinvgaussian/;
function cdfinvgaus /stolib.cdfinvgaussian/;
function icdfinvgaus /stolib.icdfinvgaussian/;
function drayl /stolib.Drayleigh /;
function pdfrayl /stolib.pdfrayleigh /;
function cdfrayl /stolib.cdfrayleigh /;
function icdfrayl /stolib.icdfrayleigh /;
function dcauchy /stolib.Dcauchy /;
function pdfcauchy /stolib.pdfcauchy /;
function cdfcauchy /stolib.cdfcauchy /;
function icdfcauchy /stolib.icdfcauchy /;
function dlognorm /stolib.Dlognormal /;
function pdflognorm /stolib.pdflognormal /;
function cdflognorm /stolib.cdflognormal /;
function icdflognorm /stolib.icdflognormal /;
function dexpo /stolib.Dexponential /;
function pdfexpo /stolib.pdfexponential/;
function cdfexpo /stolib.cdfexponential/;
function icdfexpo /stolib.icdfexponential/;
function dlogist /stolib.Dlogistic /;
function pdflogist /stolib.pdflogistic /;
function cdflogist /stolib.cdflogistic /;
function icdflogist /stolib.icdflogistic /;
function dgamma /stolib.Dgamma /;
function pdfgamma /stolib.pdfgamma /;
function cdfgamma /stolib.cdfgamma /;
function icdfgamma /stolib.icdfgamma /;
function dchisq /stolib.Dchisquare /;
function pdfchisq /stolib.pdfchisquare /;
function cdfchisq /stolib.cdfchisquare /;
function icdfchisq /stolib.icdfchisquare /;
function dweibull /stolib.Dweibull /;
function pdfweibull /stolib.pdfweibull /;
function cdfweibull /stolib.cdfweibull /;
function icdfweibull /stolib.icdfweibull /;
function dbeta /stolib.Dbeta /;
function pdfbeta /stolib.pdfbeta /;
function cdfbeta /stolib.cdfbeta /;
function icdfbeta /stolib.icdfbeta /;
function df /stolib.DF /;
function pdff /stolib.pdfF /;
function cdff /stolib.cdfF /;
function icdff /stolib.icdfF /;
function dstudt /stolib.Dstudentt /;
function pdfstudt /stolib.pdfstudentt /;
function cdfstudt /stolib.cdfstudentt /;
function icdfstudt /stolib.icdfstudentt /;
function dpareto /stolib.Dpareto /;
function pdfpareto /stolib.pdfpareto /;
function cdfpareto /stolib.cdfpareto /;
function icdfpareto /stolib.icdfpareto /;
function dgumbel /stolib.Dgumbel /;
function pdfgumbel /stolib.pdfgumbel /;
function cdfgumbel /stolib.cdfgumbel /;
function icdfgumbel /stolib.icdfgumbel /;
function dlaplace /stolib.Dlaplace /;
function pdflaplace /stolib.pdflaplace /;
function cdflaplace /stolib.cdflaplace /;
function icdflaplace /stolib.icdflaplace /;
function dtri /stolib.Dtriangular /;
function pdftri /stolib.pdftriangular /;
function cdftri /stolib.cdftriangular /;
function icdftri /stolib.icdftriangular /;
function dunifint /stolib.Duniformint /;
function cdfunifint /stolib.cdfuniformint /;
function icdfunifint /stolib.icdfuniformint /;
function pdfunifint /stolib.pdfuniformint /;
function dbino /stolib.Dbinomial /;
function cdfbino /stolib.cdfbinomial /;
function icdfbino /stolib.icdfbinomial /;
function pdfbino /stolib.pdfbinomial /;
function dnegbino /stolib.Dnegbinomial /;
function cdfnegbino /stolib.cdfnegbinomial/;
function icdfnegbino /stolib.icdfnegbinomial/;
function pdfnegbino /stolib.pdfnegbinomial/;
function dgeomet /stolib.Dgeometric /;
function cdfgeomet /stolib.cdfgeometric /;
function icdfgeomet /stolib.icdfgeometric /;
function pdfgeomet /stolib.pdfgeometric /;
function dhypergeo /stolib.Dhypergeo /;
function cdfhypergeo /stolib.cdfhypergeo /;
function icdfhypergeo /stolib.icdfhypergeo /;
function pdfhypergeo /stolib.pdfhypergeo /;
function dlogarithmic /stolib.Dlogarithmic /;
function cdflogarithmic /stolib.cdflogarithmic/;
function icdflogarithmic /stolib.icdflogarithmic/;
function pdflogarithmic /stolib.pdflogarithmic/;
function dpoisson /stolib.Dpoisson /;
function cdfpoisson /stolib.cdfpoisson /;
function icdfpoisson /stolib.icdfpoisson /;
function pdfpoisson /stolib.pdfpoisson /;
file fx; put fx;
scalar x;
*Dummy return
x = dseed(12345);
option FDDelta=1e-3;
set l /0*100/
skipg(l)
skipig(l)
k /0*100/
j /1*100/
i /1*100000/;
parameter st(l), ist(l);
skipg(l) = no;
skipig(l) = no;
st(l) = 0.1*(ord(l)-1);
ist(l) = 0.01*(ord(l)-1);
$onEchoV > cont.gms
$if %Add0% == def $set Add0 1
$if %icdfTol% == def $set icdfTol 1e-8
$if %histTol% == def $set histTol 1e-2
$if %gradTol% == def $set gradTol 1e-3
parameter histo%1%2(j)
pdf%1%2(l)
cdf%1%2(l)
cdfhisto%1%2(l)
icdf%1%2(l)
dpdf%1%2(l,*)
dcdf%1%2(l,*)
dicdf%1%2(l,*)
err0%1%2(l) 'cdf(icdf(x)) <> x'
err1%1%2(l) 'icdf(cdf(x)) <> x'
err2%1%2(l) 'cdf(x) <> histo(j<=x)'
err3%1%2(l,*) 'pdf: grad/hess <> gradn/hessn'
err4%1%2(l,*) 'cdf: grad/hess <> gradn/hessn'
err5%1%2(l,*) 'icdf: grad/hess <> gradn/hessn';
histo%1%2(j) = 0;
loop(i,
x = d%1(%3);
loop(j$sameas(j,'1'),
histo%1%2(j+round(x/st('1')-0.5))= histo%1%2(j+round(x/st('1')-0.5))+1;
)
);
pdf%1%2(l)$(not skipg(l)) = pdf%1(st(l),%3);
icdf%1%2(l)$(ist(l)>0 and ist(l)<1 and not skipig(l)) = icdf%1(ist(l),%3);
cdf%1%2(l)$(ist(l)>0 and ist(l)<1 and pdf%1%2(l)>1e-8 and not skipig(l)) = cdf%1(icdf%1%2(l),%3);
err0%1%2(l)$(ist(l)>0 and ist(l)<1 and pdf%1%2(l)>1e-8 and not skipig(l)) = cdf%1%2(l) - ist(l);
err0%1%2(l)$(abs(err0%1%2(l)) <= %icdfTol%) = 0;
cdf%1%2(l)$(not skipg(l)) = cdf%1(st(l),%3);
icdf%1%2(l)$(cdf%1%2(l)>0 and cdf%1%2(l)<1 and pdf%1%2(l)>1e-8) = icdf%1(cdf%1%2(l),%3);
err1%1%2(l)$(cdf%1%2(l)>0 and cdf%1%2(l)<1 and pdf%1%2(l)>1e-8) = icdf%1%2(l) - st(l);
err1%1%2(l)$(abs(err1%1%2(l)) <= %icdfTol%) = 0;
cdfhisto%1%2(l)$(not skipg(l)) = cdf%1(0,%3)$%Add0% + sum(j$(ord(j)*st('1')<=st(l)),histo%1%2(j))/card(i);
err2%1%2(l) = cdf%1%2(l) - cdfhisto%1%2(l);
err2%1%2(l)$(abs(err2%1%2(l)) <= %histTol%) = 0;
loop {l$(not skipg(l)),
dpdf%1%2(l,'f ') = pdf%1 ( st(l), %3);
dpdf%1%2(l,'gn') = pdf%1.gradn(1: st(l), %3);
dpdf%1%2(l,'g' ) = pdf%1.grad (1: st(l), %3);
dpdf%1%2(l,'hn') = pdf%1.hessn(1: st(l), %3);
dpdf%1%2(l,'h' ) = pdf%1.hess (1: st(l), %3);
dcdf%1%2(l,'f ') = cdf%1 ( st(l), %3);
dcdf%1%2(l,'gn') = cdf%1.gradn(1: st(l), %3);
dcdf%1%2(l,'g' ) = cdf%1.grad (1: st(l), %3);
dcdf%1%2(l,'hn') = cdf%1.hessn(1: st(l), %3);
dcdf%1%2(l,'h' ) = cdf%1.hess (1: st(l), %3);
* data(T, 'rc') = 0;
* data(T, 'ec') = 0;
};
err3%1%2(l,'g') = abs(dpdf%1%2(l,'gn') - dpdf%1%2(l,'g'))/max(1,abs(dpdf%1%2(l,'g')));
err3%1%2(l,'g')$(err3%1%2(l,'g') <= %gradTol%) = 0;
err3%1%2(l,'h') = abs(dpdf%1%2(l,'hn') - dpdf%1%2(l,'h'))/max(1,abs(dpdf%1%2(l,'h')));
err3%1%2(l,'h')$(err3%1%2(l,'h') <= %gradTol%) = 0;
err4%1%2(l,'g') = abs(dcdf%1%2(l,'gn') - dcdf%1%2(l,'g'))/max(1,abs(dcdf%1%2(l,'g')));
err4%1%2(l,'g')$(err4%1%2(l,'g') <= %gradTol%) = 0;
err4%1%2(l,'h') = abs(dcdf%1%2(l,'hn') - dcdf%1%2(l,'h'))/max(1,abs(dcdf%1%2(l,'h')));
err4%1%2(l,'h')$(err4%1%2(l,'h') <= %gradTol%) = 0;
loop {l$(ist(l)>0.05 and ist(l)<0.95 and not skipig(l)),
dicdf%1%2(l,'f ') = icdf%1 ( ist(l), %3);
dicdf%1%2(l,'gn') = icdf%1.gradn(1: ist(l), %3);
dicdf%1%2(l,'g' ) = icdf%1.grad (1: ist(l), %3);
dicdf%1%2(l,'hn') = icdf%1.hessn(1: ist(l), %3);
dicdf%1%2(l,'h' ) = icdf%1.hess (1: ist(l), %3);
};
err5%1%2(l,'g') = abs(dicdf%1%2(l,'gn') - dicdf%1%2(l,'g'))/max(1,abs(dicdf%1%2(l,'g')));
err5%1%2(l,'g')$(err5%1%2(l,'g') <= %gradTol%) = 0;
err5%1%2(l,'h') = abs(dicdf%1%2(l,'hn') - dicdf%1%2(l,'h'))/max(1,abs(dicdf%1%2(l,'h')));
err5%1%2(l,'h')$(err5%1%2(l,'h') <= %gradTol%) = 0;
abort$(card(err0%1%2) + card(err1%1%2) + card(err2%1%2) + card(err3%1%2) + card(err4%1%2) + card(err5%1%2))
err0%1%2, err1%1%2, err2%1%2, err3%1%2, err4%1%2, err5%1%2;
$offEcho
$onEchoV > disc.gms
$if %icdfTol% == def $set icdfTol 0
$if %histTol% == def $set histTol 1e-2
parameter histo%1%2(k)
pdf%1%2(l)
cdf%1%2(l)
cdfhisto%1%2(l)
icdf%1%2(l)
err0%1%2(l) 'cdf(icdf(x)) <> x'
err1%1%2(l) 'icdf(cdf(x)) <> x'
err2%1%2(l) 'cdf(x) <> histo(j<=x)';
histo%1%2(k) = 0;
loop(i,
x = d%1(%3);
loop(k$sameas(k,'0'),
histo%1%2(k+x)= histo%1%2(k+x)+1;
)
);
pdf%1%2(l)$(not skipg(l)) = pdf%1(st(l),%3);
icdf%1%2(l)$(ist(l)>0 and ist(l)<1 and not skipig(l)) = icdf%1(ist(l),%3);
cdf%1%2(l)$(ist(l)>0 and ist(l)<1 and pdf%1%2(l)>1e-8 and not skipig(l)) = cdf%1(icdf%1%2(l),%3);
err0%1%2(l)$(ist(l)>0 and ist(l)<1 and pdf%1%2(l)>1e-8 and not skipig(l)) = icdf%1(cdf%1%2(l),%3) - icdf%1(ist(l),%3);
err0%1%2(l)$(abs(err0%1%2(l)) <= %icdfTol%) = 0;
cdf%1%2(l)$(not skipg(l)) = cdf%1(st(l),%3);
icdf%1%2(l)$(cdf%1%2(l)>0 and cdf%1%2(l)<1 and pdf%1%2(l)>1e-8) = icdf%1(cdf%1%2(l),%3);
err1%1%2(l)$(cdf%1%2(l)>0 and cdf%1%2(l)<1 and pdf%1%2(l)>1e-8) = icdf%1%2(l) - st(l);
err1%1%2(l)$(abs(err1%1%2(l)) <= %icdfTol%) = 0;
cdfhisto%1%2(l)$(not skipg(l)) = sum(k$(ord(k)-1<=st(l)),histo%1%2(k))/card(i);
err2%1%2(l) = cdf%1%2(l) - cdfhisto%1%2(l);
err2%1%2(l)$(abs(err2%1%2(l)) <= %histTol%) = 0;
abort$(card(err0%1%2) + card(err1%1%2) + card(err2%1%2))
err0%1%2, err1%1%2, err2%1%2;
$offEcho
$set Add0 def
$set icdfTol def
$set histTol def
$set gradTol def
st(l) = 1*(ord(l)-1);
put_utility 'log' / 'Uniform Integer';
$batInclude disc unifint _1_5 1,5
$batInclude disc unifint _3_15 3,15
put_utility 'log' / 'Binomial';
$batInclude disc bino _20_05 20,0.5
$batInclude disc bino _20_07 20,0.7
$batInclude disc bino _40_05 40,0.5
put_utility 'log' / 'Negative Binomial';
$batInclude disc negbino _10_02 10,0.2
$batInclude disc negbino _10_05 10,0.5
$batInclude disc negbino _10_08 10,0.8
put_utility 'log' / 'Geometric';
$batInclude disc geomet _02 0.2
$batInclude disc geomet _05 0.5
$batInclude disc geomet _08 0.8
put_utility 'log' / 'Hypergeometric';
$batInclude disc hypergeo _30_20_20 30,20,20
$batInclude disc hypergeo _60_50_20 60,50,20
$batInclude disc hypergeo _60_20_20 60,20,20
skipg(l)$(st(l)<1)=yes;
put_utility 'log' / 'Logarithmic';
$batInclude disc logarithmic _033 0.33
$batInclude disc logarithmic _066 0.66
$batInclude disc logarithmic _099 0.99
skipg(l)=no;
put_utility 'log' / 'Poisson';
$batInclude disc poisson _1 1
$batInclude disc poisson _4 4
$batInclude disc poisson _10 10
st(l) = 0.1*(ord(l)-1);
skipg(l)$(st(l)=1)=yes;
skipg(l)$(st(l)=5)=yes;
put_utility 'log' / 'Uniform';
$batInclude cont unif _1_5 1,5
skipg(l)=no;
put_utility 'log' / 'Normal';
$batInclude cont norm _0 0,1
$batInclude cont norm _5 5,1
skipg(l)$(st(l)<=0.1)=yes;
$set Add0 0
$set icdfTol 1e-4
put_utility 'log' / 'Inverse Gaussian';
$batInclude cont invgaus _2_1 2,1
$batInclude cont invgaus _2_5 2,5
$batInclude cont invgaus _2_30 2,30
$batInclude cont invgaus _1_10 1,1
$batInclude cont invgaus _1_02 1,0.2
$batInclude cont invgaus _1_30 1,3
$batInclude cont invgaus _3_10 3,1
$batInclude cont invgaus _3_02 3,0.2
skipg(l)=no;
$set Add0 def
$set icdfTol def
skipg(l)$(st(l)=0)=yes;
put_utility 'log' / 'Rayleigh';
$batInclude cont rayl _05 0.5
$batInclude cont rayl _10 1
$batInclude cont rayl _20 2
$batInclude cont rayl _30 3
$batInclude cont rayl _40 4
skipg(l)=no;
skipig(l)$(ist(l)<=0.06 or ist(l)>=0.94)=yes;
put_utility 'log' / 'Cauchy';
$batInclude cont cauchy _0_05 0,0.5
$batInclude cont cauchy _0_10 0,1.0
$batInclude cont cauchy _0_20 0,2.0
$batInclude cont cauchy _2_10 2,1.0
$batInclude cont cauchy _5_10 5,1.0
skipig(l)=no;
$set Add0 0
$set gradTol 1e-2
skipg(l)$(st(l)<=0)=yes;
put_utility 'log' / 'Lognormal';
$batInclude cont lognorm _0_025 0,0.25
$batInclude cont lognorm _0_1 0,1
$batInclude cont lognorm _0_10 0,10
skipg(l)=no;
$set Add0 def
$set gradTol def
skipg(l)$(st(l)=0)=yes;
put_utility 'log' / 'Exponential';
$batInclude cont expo _05 0.5
$batInclude cont expo _10 1.0
$batInclude cont expo _15 1.5
skipg(l)=no;
put_utility 'log' / 'Logistic';
$batInclude cont logist _5_2 5,2
$batInclude cont logist _9_3 9,3
$batInclude cont logist _9_4 9,4
$batInclude cont logist _6_2 6,2
$batInclude cont logist _2_1 2,1
skipg(l)$(st(l)=0)=yes;
put_utility 'log' / 'Gamma';
$batInclude cont gamma _1_20 1,2
$batInclude cont gamma _2_20 2,2
$batInclude cont gamma _3_20 3,2
$batInclude cont gamma _5_10 5,1
$batInclude cont gamma _9_05 9,0.5
skipg(l)=no;
skipg(l)$(st(l)=0)=yes;
put_utility 'log' / 'ChiSquare';
$batInclude cont chisq _1 1
$batInclude cont chisq _2 2
$batInclude cont chisq _3 3
skipg(l)=no;
skipg(l)$(st(l)=0)=yes;
put_utility 'log' / 'Weibull';
$batInclude cont weibull _05_1 0.5,1
$batInclude cont weibull _10_1 1,1
$batInclude cont weibull _15_1 1.5,1
$batInclude cont weibull _50_1 5,1
$batInclude cont weibull _50_2 5,2
skipg(l)=no;
st(l) = 0.01*(ord(l)-1);
skipg(l)$(st(l)<=0.1 or st(l)>=0.9)=yes;
put_utility 'log' / 'Beta';
$batInclude cont beta _05_05 0.5,0.5
$batInclude cont beta _50_10 5,1
$batInclude cont beta _10_30 1,3
$batInclude cont beta _20_20 2,2
$batInclude cont beta _20_50 2,5
st(l) = 0.1*(ord(l)-1);
skipg(l)=no;
skipg(l)$(st(l)=0)=yes;
skipig(l)$(ist(l)<=0.1 or ist(l)>=0.9)=yes;
put_utility 'log' / 'F';
$batInclude cont f _1_1 1,1
$batInclude cont f _2_1 2,1
$batInclude cont f _5_2 5,2
$batInclude cont f _100_1 100,1
$batInclude cont f _100_100 100,100
skipg(l)=no;
skipig(l)=no;
skipig(l)$(ist(l)<=0.06 or ist(l)>=0.94)=yes;
put_utility 'log' / 'Student T';
$batInclude cont studt _1 1
$batInclude cont studt _2 2
$batInclude cont studt _5 5
skipig(l)=no;
$set Add0 0
skipg(l)$(st(l)<=1)=yes;
skipig(l)$(ist(l)<=0.06 or ist(l)>=0.94)=yes;
put_utility 'log' / 'Pareto';
$batInclude cont pareto _1_1 1,1 $(st(l)>=1)
$batInclude cont pareto _1_2 1,2 $(st(l)>=1)
$batInclude cont pareto _1_3 1,3 $(st(l)>=1)
skipg(l)=no;
skipg(l)$(st(l)<=2)=yes;
$batInclude cont pareto _2_2 2,2 $(st(l)>=2)
skipg(l)=no;
skipig(l)=no;
$set Add0 def
put_utility 'log' / 'Gumbel';
$batInclude cont gumbel _05_2 0.5,2
$batInclude cont gumbel _10_2 1,2
$batInclude cont gumbel _15_3 1.5,3
$batInclude cont gumbel _30_4 3,4
skipg (l)$( st(l)=0 )=yes;
skipig(l)$(ist(l)=0.5)=yes;
put_utility 'log' / 'Laplace';
$batInclude cont laplace _0_1 0,1
$batInclude cont laplace _0_2 0,2
$batInclude cont laplace _0_4 0,4
skipg(l)=no;
skipg(l)$(st(l)=6)=yes;
$batInclude cont laplace _6_2 6,2
skipg(l)=no;
skipg(l)$(st(l)=5)=yes;
$batInclude cont laplace _5_4 5,4
skipg (l)=no;
skipig(l)=no;
skipg(l)$(st(l)=1)=yes;
skipg(l)$(st(l)=2)=yes;
skipg(l)$(st(l)=5)=yes;
skipig(l)$(ist(l)=(2-1)/(5-1))=yes;
put_utility 'log' / 'Triangular';
$batInclude cont tri _1_2_5 1,2,5
skipg(l)=no;
skipig(l)=no;
skipg(l)$(st(l)=3)=yes;
skipg(l)$(st(l)=6)=yes;
skipg(l)$(st(l)=9)=yes;
skipig(l)$(ist(l)=(6-3)/(9-3))=yes;
$batInclude cont tri _3_6_9 3,6,9
skipg(l)=no;
skipig(l)=no;
skipg(l)$(st(l)=2)=yes;
skipg(l)$(st(l)=7)=yes;
$batInclude cont tri _2_7_7 2,7,7
skipg(l)=no;
skipg(l)$(st(l)=8)=yes;
skipg(l)$(st(l)=10)=yes;
$batInclude cont tri _8_8_10 8,8,10
skipg(l)=no;
$call rm -rf disc.gms cont.gms