1:- module(rint, [pbinom/5]).    2
    3:- multifile r_hook/1.    4
    5:- reexport(interval).    6:- reexport(r), r_initialize.

   14% Binomial distribution
   15% - Cumulated density lower-tail: interval(pbinom(X, N, P, true), Res)
   16% - Cumulated density upper-tail: interval(pbinom(X, N, P, false), Res)
   17% - Quantile: interval(qbinom(Alpha, N, P, true), Res)
   18% - Density: interval(dbinom(X, N, P), Res)
   19%
   20% Normal distribution
   21% - Cumulated density: interval(pnorm(X, Mu, Sigma), Res)
   22% - Quantile: interval(qnorm(P, Mu, Sigma), Res)
   23% - Density: interval(dnorm(X, Mu, Sigma), Res)
   24
   25%
   26% Skip R vectors
   27%
   28interval:int_hook(:, colon(_, _)).
   29colon(A, A).
   30
   31%
   32% Obtain atoms or functions from R
   33%
   34interval:eval_hook(Atom, Res) :-
   35    atomic(Atom),
   36    r_hook(Atom),
   37    !,
   38    r(Atom, Res).
   39
   40interval:eval_hook(Expr, Res) :-
   41    compound(Expr),
   42    compound_name_arity(Expr, Name, Arity),
   43    r_hook(Name/Arity),
   44    !,
   45    r(Expr, Res).
   46
   47r_hook(true).
   48r_hook(false).
   49
   50%
   51% Binomial distribution
   52%
   53interval:int_hook(pbinom, pbinom(..., ..., ..., atomic)).
   54
   55% lower tail
   56pbinom(X, N, P, atomic(true), Res) :-
   57    !,
   58    interval(pbinom0(X, N, P), Res).
   59
   60% upper tail
   61pbinom(X, N, P, atomic(false), Res) :-
   62    interval(pbinom1(X, N, P), Res).
   63
   64r_hook(pbinom0/3).
   65interval:mono(pbinom0/3, [+, -, -]).
   66
   67r_hook(pbinom1/3).
   68interval:mono(pbinom1/3, [-, +, +]).
   69
   70%
   71% Quantile function - hier weiter
   72%
   73interval:int_hook(qbinom/4, []).
   74
   75% lower tail
   76interval:int_hook(qbinom(Alpha, N, P, true), Res, Opt) :-
   77    !,
   78    interval(qbinom0(Alpha, N, P), Res, Opt).
   79
   80r_hook(qbinom0/3).
   81interval:mono(qbinom0/3, [+, +, +]).
   82
   83% upper tail
   84interval:int_hook(qbinom(Alpha, N, P, false), Res, Opt) :-
   85    interval(qbinom1(Alpha, N, P), Res, Opt).
   86
   87r_hook(qbinom1/3).
   88interval:mono(qbinom1/3, [-, +, +]).
   89
   90%
   91% Density
   92%
   93interval:int_hook(dbinom/3, []).
   94
   95% left to X / N
   96interval:int_hook(dbinom(X1...X2, N1...N2, P1...P2), Res, Opt) :-
   97    X2 < N1 * P1,
   98    !,
   99    interval(dbinom0(X1...X2, N1...N2, P1...P2), Res, Opt).
  100
  101r_hook(dbinom0/3).
  102interval:mono(dbinom0/3, [+, -, -]).
  103
  104% right to X / N
  105interval:int_hook(dbinom(X1...X2, N1...N2, P1...P2), Res, Opt) :-
  106    X1 > N2 * P2,
  107    !,
  108    interval(dbinom1(X1...X2, N1...N2, P1...P2), Res, Opt).
  109
  110r_hook(dbinom1/3).
  111interval:mono(dbinom1/3, [-, +, +]).
  112
  113% otherwise
  114interval:int_hook(dbinom(X1...X2, N1...N2, P1...P2), Res, _) :-
  115    r(dbinom2(X1, X2, N1, N2, P1, P2), #(L, U)),
  116    Res = L...U.
  117
  118%
  119% Normal distribution
  120%
  121r_hook(pnorm0/1).
  122interval:mono(pnorm0/1, [+]).
  123
  124interval:int_hook(pnorm/3, []).
  125interval:int_hook(pnorm(X, Mu, Sigma), Res, Opt) :-
  126     interval((X - Mu)/Sigma, Z, Opt),
  127     interval(pnorm0(Z), Res, Opt).
  128
  129%
  130% Quantile function
  131%
  132r_hook(qnorm0/1).
  133interval:mono(qnorm0/1, [+]).
  134
  135interval:int_hook(qnorm/3, []).
  136interval:int_hook(qnorm(P, Mu, Sigma), Res, Opt) :-
  137     interval(qnorm0(P), Z, Opt),
  138     interval(Mu + Z * Sigma, Res, Opt).
  139
  140%
  141% Density
  142%
  143r_hook(dnorm1/1).
  144interval:mono(dnorm1/1, [+]).
  145
  146r_hook(dnorm2/1).
  147interval:mono(dnorm2/1, [-]).
  148
  149interval:int_hook(dnorm/3, []).
  150interval:int_hook(dnorm(X, Mu, Sigma), Res, Opt) :-
  151    interval((X - Mu)/Sigma, Z, Opt),
  152    interval(1/Sigma * dnorm0(Z), Res, Opt).
  153
  154interval:int_hook(dnorm0/1, []).
  155interval:int_hook(dnorm0(A...B), Res, Opt) :-
  156    B =< 0,
  157    !,
  158    interval(dnorm1(A...B), Res, Opt).
  159
  160interval:int_hook(dnorm0(A...B), Res, Opt) :-
  161    A >= 0,
  162    !,
  163    interval(dnorm2(A...B), Res, Opt).
  164
  165% mixed
  166interval:int_hook(dnorm0(A