mercury-rmath: Mercury library for the standalone Rmath library

Introduction

This is an alpha quality Mercury library for using the standalone Rmath library. The library is based on C and the installation is currently only for the ast_fast.gc grade.

This has been tested on Linux with the Rmath library installed. On Ubuntu/Debian, you would need apt install r-mathlib. For the lib and cflags paths below, I have also used pkg-config.

Following the R language, I have used a GPL licence.

Feedback is welcome.

Installation

For an installation in the local directory, use make install, which is currently defined as:

rm -rf Mercury lib # is this required?
mmc --make  --c-include-directory `pkg-config --cflags libRmath` `pkg-config --libs libRmath` --no-libgrade --libgrade asm_fast.gc --install-prefix . librmath.install

To install to the main Mercury library, remove: --install-prefix .

Documentation

See the interface in https://github.com/mclements/mercury-rmath/blob/main/rmath.m.

Test

As a test file:

cat test_rmath.m

:- module test_rmath.

:- interface.
:- import_module io.
:- pred main(io::di, io::uo) is det.

:- implementation.
:- import_module int, float, rmath, pair, bool.

:- type alternative ---> two_sided ; less ; greater.
:- pred poisson_ci(float::in, float::in, alternative::in, pair(float)::out).
poisson_ci(X, Conflevel, Alternative, Interval) :-
    Alpha = (1.0-Conflevel)/2.0,
    Pl = (func(Xi,Alphai) = (if Xi=0.0 then 0.0 else rmath.qgamma(Alphai,Xi, 1.0, 1, 0))),
    Pu = (func(Xi,Alphai) = rmath.qgamma(1.0-Alphai, Xi+1.0, 1.0, 1, 0)),
    Interval = (Alternative = less -> (0.0 - Pu(X, 1.0-Conflevel))
	       ;
	       Alternative = greater -> (Pl(X, 1.0-Conflevel) - 1.0)
	       ;
	       %% two_sided
	       (Pl(X,Alpha) - Pu(X, Alpha))).

:- func for_loop(func(int,int) = int, int, int, int) = int.
for_loop(Fun, I, Finish, Agg) = (if I>Finish then Agg else for_loop(Fun, I+1, Finish, Fun(I,Agg))).
:- func count(func(int) = bool, int, int) = int.
count(Predicate, Start, Finish) = Result :-
    Result = for_loop(func(I, Y) = (if Predicate(I)=yes then Y+1 else Y), Start, Finish, 0).

:- func loop1(int,float,float) = int.
loop1(Ni,M,D) = (if rmath.dpois(float(Ni),M,0)>D then loop1(Ni*2,M,D) else Ni).
:- func poisson_test(float,float,float,alternative) = float.
poisson_test(X, T, R, Alternative) = Result :-
    M = R*T,
    (Alternative = less -> Result = rmath.ppois(X,M,1,0)
    ;
    Alternative = greater -> Result = rmath.ppois(X-1.0,M,0,0)
    ;
    %% Alternative = two_sided
    (M = 0.0 -> Result = (X=0.0 -> 1.0; 0.0)
	   ;
	   (Relerr = 1.00000001,
            D = rmath.dpois(X,M,0),
	    Dstar = D * Relerr,
	    Pred = (func(I) = (if rmath.dpois(float(I),M,0) =< Dstar then yes else no)),
            (X=M -> Result = 1.0
	     ;
	     X<M ->
	     (N = loop1(ceiling_to_int(2.0*M-X),M,D),
	      Y = count(Pred, ceiling_to_int(M), N),
	      Result = rmath.ppois(X,M,1,0) + rmath.ppois(float(N)-float(Y),M,0,0))
	     ;
	     %% X>M
	     (Y = count(Pred,0,floor_to_int(M)),
	      Result = rmath.ppois(float(Y)-1.0,M,1,0) + rmath.ppois(X-1.0, M,0,0)))))).

main(!IO) :-
    set_seed(1219005190u32, 1225564623u32, !IO),
    runif(0.0, 1.0, U, !IO),
    runif(0.0, 1.0, U2, !IO),
    poisson_ci(5.0, 0.95, two_sided, Interval),
    P = poisson_test(5.0, 1.0, 1.0, two_sided),
    P2 = poisson_test(5.0, 10.0, 1.0, two_sided),
    io.write_line({U, U2,
		   m_e,
		   rmath.pnorm(1.96, 0.0, 1.0, 1, 0),
		   rmath.qnorm(0.975, 0.0, 1.0, 1, 0),
		   Interval,
		   P, P2
		  },
		  !IO).

This can be run using make test, which is currently defined as:

mmc --make --mld ./lib/mercury --ml rmath `pkg-config --libs libRmath` test_rmath && ./test_rmath

{0.6801521621365453, 0.5965995308003852, 2.718281828459045, 0.9750021048517796, 1.9599639845400536, 1.6234863901184209 - 11.668332079322667, 0.0036598468273437144, 0.15054443581369478}

To compare with R output:

RNGkind("Marsaglia-Multicarry")
set.seed(3)
## .Random.seed = c(10401L, 1219005190L, 1225564623L)
P = poisson.test(5)
P2 = poisson.test(5,10)
c(runif(2), exp(1), pnorm(1.96), qnorm(0.975), P$conf.int[1], P$conf.int[2], P$p.value, P2$p.value)

0.680152162136545
0.596599530800385
2.71828182845905
0.97500210485178
1.95996398454005
1.62348639011842
11.6683320793227
0.00365984682734371
0.150544435813695