rimパッケージでRからMaximaを使う – 伊東生態統計研究室

RからMaximaを使えるようにするrimというパッケージがあるということでためしてみました。

準備

rimパッケージはCRANからインストールできます。

install.packages("rim")

使うにはlibrary関数で呼び出します。フォーマットにlatexを指定します。

library("rim")
## Maxima successfully registered as knitr engine!

maxima.options(engine.format = "latex", 
               engine.label = TRUE,
               inline.format = "latex", 
               inline.label = FALSE)

使用法

RStudioやQuartoでは、チャンクの先頭行を```{maxima}とすれば、Maximaのコードとして解釈されます。

Maximaで正規分布の密度関数を定義してみました。

(%i1) f(x,mu,sigma) := 1/sqrt(2*%pi*sigma^2)*
  exp(-(x-mu)^2/(2*sigma^2));

\[\mathtt{(\textit{%o}_{1})}\quad f\left(x , \mu , \sigma\right):=\frac{1}{\sqrt{2\,\pi\,\sigma^2}}\,\exp \left(\frac{-\left(x-\mu\right)^2}{2\,\sigma^2}\right)\]

ラベルに(%o1)と表示されるはずですが、%oがコメントと解釈されてしまったようです。これはとりあえず目をつむっておきます。

次に、平均0、標準偏差1の正規分布でx=1のときの密度です。

(%i2) float(f(1, 0, 1));

\[\mathtt{(\textit{%o}_{2})}\quad 0.24197072451914334\]

Rで計算した値と比較してみます。

dnorm(1, 0, 1)
## [1] 0.2419707

\(\int_{-\infty}^{ 0}N(x \mid 0, 1)dx\)を求めてみます。

(%i3) integrate(f(x, 0, 1), x, minf, 0);

\[\mathtt{(\textit{%o}_{3})}\quad \frac{1}{2}\]

Rコンソールで使う

RStudioではMaximaの出力が表示されますが、Rコンソールでは何も表示されません。

m1 <- maxima.get("diff((x^2 + 1) * log(x), x)")
print(m1)
## (%o1) 2*x*log(x)+(x^2+1)/x

この記事はQuartoでレンダリングしているので、m1の値が表示されていますが、コードをコピーしてコンソールで実行しても何も表示されません。

print(m1)

str関数で構造を確認します。

str(m1)
## List of 2
##  $ wtl:List of 5
##   ..$ linear: chr "(%o1) 2*x*log(x)+(x^2+1)/x"
##   ..$ ascii : chr [1:4] "                                            2" "                                           x  + 1" "(%o1)                         2 x log(x) + ------" "                                             x"
##   ..$ latex : chr "$$\\mathtt{(\\textit{\\%o}_{1})}\\quad 2\\,x\\,\\log x+\\frac{x^2+1}{x}$$"
##   ..$ inline: chr "$\\mathtt{(\\textit{\\%o}_{1})}\\quad 2\\,x\\,\\log x+\\frac{x^2+1}{x}$"
##   ..$ mathml: chr [1:7] "<math xmlns=\"http://www.w3.org/1998/Math/MathML\"> <mi>mlabel</mi> " " <mfenced separators=\"\"><msub><mi>%o</mi> <mn>1</mn></msub> " " <mo>,</mo><mn>2</mn> <mspace width=\"thinmathspace\"/><mi>x</mi> " " <mspace width=\"thinmathspace\"/><mi>log</mi> <mi>x</mi> <mo>+</mo> " ...
##  $ wol:List of 6
##   ..$ linear: chr "2*x*log(x)+(x^2+1)/x"
##   ..$ ascii : chr [1:4] "              2" "             x  + 1" "2 x log(x) + ------" "               x"
##   ..$ latex : chr "$$2\\,x\\,\\log x+\\frac{x^2+1}{x}$$"
##   ..$ inline: chr "$2\\,x\\,\\log x+\\frac{x^2+1}{x}$"
##   ..$ mathml: chr [1:5] "<math xmlns=\"http://www.w3.org/1998/Math/MathML\"> " " <mn>2</mn> <mspace width=\"thinmathspace\"/><mi>x</mi> " " <mspace width=\"thinmathspace\"/><mi>log</mi> <mi>x</mi> <mo>+</mo> " " <mfrac><mrow><msup><mrow><mi>x</mi> </mrow> <mn>2</mn> </msup> " ...
##   ..$ rstr  : chr "(((x ^ -1L) * (1L + (x ^ 2L))) + (2L * x * log(x)))"
##  - attr(*, "input.label")= chr "%i1"
##  - attr(*, "output.label")= chr "%o1"
##  - attr(*, "command")= chr "diff((x^2 + 1) * log(x), x)"
##  - attr(*, "suppressed")= logi FALSE
##  - attr(*, "parsed")= language, mode "(": (((x^-1L) * (1L + (x^2L))) + (2L * x * log(x)))
##  - attr(*, "from_engine")= logi FALSE
##  - attr(*, "class")= chr "maxima"

wtlがknitrエンジン用の出力、wolがインライン用の出力のようです。

以下のようにすると、コンソールでもインライン形式でのlinear出力を表示できます。

print(m1$wol$linear)
## [1] "2*x*log(x)+(x^2+1)/x"

以下のようにすると、同じくアスキーが表示されます。

cat(m1$wol$ascii, sep = "\n")
##               2
##              x  + 1
## 2 x log(x) + ------
##                x

LaTeXのコードでの出力です。

m1$wol$latex
## [1] "$$2\\,x\\,\\log x+\\frac{x^2+1}{x}$$"

これは、latex2expパッケージのTeX関数でレンダリングできます。

library(latex2exp)

plot(TeX(m1$wol$latex))