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.
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 .
See the interface in https://github.com/mclements/mercury-rmath/blob/main/rmath.m.
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 |