Accessing R from maxima, WAS: Accessing maxima from the R statistical program
Subject: Accessing R from maxima, WAS: Accessing maxima from the R statistical program
From: Mario Rodriguez
Date: Sat, 25 Nov 2006 11:36:42 +0100
Dear all,
Trying to call R from Maxima, I wrote the attached toy-package.
It works as explained here:
1.- Maxima function 'Rtest_means_difference' writes a file
called 'maxout.r'
2.- R is called by executing system command
R -q --vanilla < maxout.r > trash
(the problem here is that the 'trash' file stores useless
information, but if I don't write this I get the R
session merged into the Maxima session)
3.- R writes file 'rout.maxima' which is a column with
the results computed by R.
4.- 'Rtest_means_difference' calls 'read_list' of the
numericalio package to read 'rout.maxima'
5.- 'Rtest_means_difference' constructs the
'inference_result' object and sends it to the console.
This procedure mimics the way Maxima calls gnuplot for plotting.
A session follows (Linux + cvs Maxima):
/********** begin Maxima session ***********/
(%i99) load("Rstats.mac")$
(%i100) Rtest_means_difference(x,y,'mu=3);
| Welch Two Sample t-test
|
| statistic = 1.56664539806187
|
| parameter = 8.62758740184604
|
| p_value = 0.153090583808592
|
(%o100) | conf_int = [- 4.85473914070372, 45.4947391407037]
|
| conf_level = 0.95
|
| estimate = [37.2, 16.88]
|
| null_value = 3
|
| alternative = twosided
(%i101) Rtest_means_difference(x,y,alternative='less,
'conflevel=99/100,'varequal='true);
| Two Sample t-test
|
| statistic = 1.76599612451501
|
| parameter = 9
|
| p_value = 0.944396790074707
|
(%o101) | conf_int = [minf, 52.7841814558168]
|
| conf_level = 0.99
|
| estimate = [37.2, 16.88]
|
| null_value = 0
|
| alternative = less
/********** end Maxima session ***********/
At least, it works.
Best wishes.
--
Mario Rodriguez Riotorto
www.biomates.net
-------------- next part --------------
/* COPYRIGHT NOTICE
Copyright (C) 2006 Mario Rodriguez Riotorto
This program is free software; you can redistribute
it and/or modify it under the terms of the
GNU General Public License as published by
the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it
will be useful, but WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details at
http://www.gnu.org/copyleft/gpl.html
*/
/* INTRODUCTION
Trying to make Maxima and R good friends ...
mario AATT edu ADOT xunta ADOT es
*/
/* LOADING PACKAGES */
load(numericalio)$
load(inference_result)$
/* M2R CONVERTERS */
/* Converts a Maxima list of numbers into an R */
/* list in string format, that is: */
/* [1,2,3] => "c(1,2,3)" */
M2Rlist(l):=block([n,str],
if not listp(l)
then return(-1),
l: map('float,l),
if not every('identity,map('numberp,l))
then return(-2),
n: length(l),
if n=0
then return("c()"),
if n=1
then return(sconcat("c(",l[1],")")),
str: sconcat("c(",l[1]),
for k:2 thru n do
str: sconcat(str,", ", l[k]),
sconcat(str,")") )$
/* Converts Maxima alternative options into */
/* R alternative options */
M2Ralternative(alt):=
if alt='twosided
then "\"two.sided\""
else if alt='greater
then "\"greater\""
else "\"less\"" $
/* Converts Maxima boolean values into */
/* R boolean values */
M2Rboolean(boo):=
if boo='true
then "TRUE"
else "FALSE" $
/* R2M CONVERTERS */
/* Converts R alternative options into */
/* Maxima alternative options */
R2Malternative(alt):=
if alt=" two.sided "
then 'twosided
else if alt=" greater "
then 'greater
else 'less $
/* Converts R infinities into Maxima infinities */
R2Minfinities(val):=
if val = - 'Inf
then 'minf
else if val = 'Inf
then 'inf
else val $
/* TEST FUNCTIONS */
/* Unpaired difference of means t-test. See 't.test' in R manual. */
/* 50% of this code is to control input arguments and options, 30% */
/* is to write the script and to call R, and the last 20% is to */
/* construct the Maxima object 'inference_result'. */
/* Admits the following options: */
/* 'alternative='twosided: this is the alternative hypothesis H1; */
/* valid values are: 'twosided, 'greater (m1>m2) and */
/* 'less (m1<m2). */
/* 'mu=0: this is the difference value to be tested */
/* 'varequal=false: whether variances are equal or not */
/* 'conflevel=0.95: confidence level, a value in (0,1) */
/* */
/* The output of this function is an 'inference_result' object */
/* with the following results: */
/* 1. 'statistic: statistic used in the procedure */
/* 2. 'parameter: degrees of freedom */
/* 3. 'p_value: p-value of the sample statistic */
/* 4. 'conf_interval= confidence interval */
/* 5. 'conf_level= confidence level */
/* 6. 'estimate= means estimates (m1, m2) */
/* 7. 'null_value: difference to be checked */
/* 8. 'alternative: alternative hypothesis */
/* Examples:
x: [20.4,62.5,61.3,44.2,11.1,23.7]$
y: [1.2,6.9,38.7,20.4,17.2]$
Rtest_means_difference(x,y,'varequal='true);
Rtest_means_difference(x,y,alternative='greater);
Rtest_means_difference(x,y,alternative='less,'conflevel=99/100);
Rtest_means_difference(x,y,'mu=3);
*/
Rtest_means_difference(x,y,[select]):=
block([aux,nx,ny,rout,options,defaults],
/* controlling arguments x and y */
if not listp(x) or not every('identity,map('numberp,x))
then error("Argument 'x must be a list of numbers"),
x: map('float,x),
nx: length(x),
if nx<2
then error("List 'x' must have more than one number"),
if not listp(y) or not every('identity,map('numberp,y))
then error("Argument 'y must be a list of numbers"),
y: map('float,y),
ny: length(y),
if ny<2
then error("List 'y' must have more than one number"),
/* updating and controlling options */
options: ['alternative, 'mu, 'varequal, 'conflevel],
defaults: ['twosided, 0, 'false, 0.95],
for i in select do(
aux: ?position(lhs(i),options),
if numberp(aux) and aux <= length(options) and aux >= 1
then defaults[aux]: rhs(i)),
if not member(defaults[1],['twosided, 'greater, 'less])
then error("Option 'alternative is not correct"),
if not numberp(defaults[2])
then error("Option 'mu' must be a number"),
defaults[2]: float(defaults[2]),
if not member(defaults[3], ['true, 'false])
then error("Option 'varequal must be a 'true or 'false"),
if not numberp(defaults[4])
then error("Option 'conflevel must be a number"),
defaults[4]: float(defaults[4]),
if defaults[4] <= 0.0 or defaults[4] >= 1.0
then error("Option 'conflevel' can't be outside interval (0,1)"),
/* Let's call R to make the job */
with_stdout("maxout.r",
print("x <- ", M2Rlist(x)),
print("y <- ", M2Rlist(y)),
print("z <- t.test(x,y,",
"alternative=", M2Ralternative(defaults[1]),",",
"mu=",defaults[2],",",
"paired= FALSE,",
"var.equal=", M2Rboolean(defaults[3]),",",
"conf.level=", defaults[4],")"),
print("output <- c(z$statistic[[1]],"),
print(" z$parameter[[1]],"),
print(" z$p.value[[1]],"),
print(" z$conf.int[[1]],"),
print(" z$conf.int[[2]],"),
print(" attributes(z$conf.int)[[1]],"),
print(" z$estimate[[1]],"),
print(" z$estimate[[2]],"),
print(" z$null.value[[1]],"),
print(" paste(\"\\\"\",z$alternative[[1]],\"\\\"\"),"),
print(" paste(\"\\\"\",z$method[[1]],\"\\\"\") )"),
print("write(output, file = \"rout.maxima\", ncolumns = 1)") ),
system("R -q --vanilla < maxout.r > trash"),
rout: read_list("rout.maxima"),
/* Let's translate some R non numeric results into Maxima */
rout[10]: R2Malternative(rout[10]),
rout[4]: R2Minfinities(rout[4]),
rout[5]: R2Minfinities(rout[5]),
/* Construct the Maxima inference_result object */
inference_result(rout[11], /* this is the method */
[ ['statistic, rout[1]],
['parameter, rout[2]],
['p_value, rout[3]],
['conf_int, [rout[4], rout[5]]],
['conf_level, rout[6]],
['estimate, [rout[7], rout[8]]],
['null_value, rout[9]],
['alternative, rout[10]] ],
[1,2,3,4,5,6,7,8]) )$