裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

RStudio での日本語の問題

2018年12月17日 | ブログラミング

Version 1.1.463 になって,解消されたようだ

コメント

RStudio で日本語が使えなくなっていた

2018年09月30日 | ブログラミング

RStudio Version 1.1.456 を立ち上げると,

During startup - Warning messages:
1: Setting LC_CTYPE failed, using "C"
2: Setting LC_COLLATE failed, using "C"
3: Setting LC_TIME failed, using "C"
4: Setting LC_MESSAGES failed, using "C"
5: Setting LC_MONETARY failed, using "C"

と出て,案の定,日本語を含むプログラムを読み込むと,文字化けしている。

もう,これは他のことで経験済みなので,慌てることなく

Sys.setlocale("LC_ALL", "ja_JP.UTF-8")

とやって,問題なし。

コメント

pythontex は,使えないな!

2018年09月04日 | RとLaTeX

pythontex は,望み薄かな?

実装が古いのと,日本語化の関連で,ちょっと面倒くさいというか無理っぽい。

 

やっぱ,R と Sweave の組み合わせが最強。

コメント

Pweave の入力と出力の区別

2018年09月04日 | RとLaTeX

そのまんまでは,入力と出力のフォーマットに全く違いがない。ので,まぎらわしい。

Sweave のようには行かず,インタラクティブモードのように,入力の前にプロンプトがつくようにはできないようだ。

そこで,LaTeX の listings を使う。

listings.sty と jlisting.sty を用意して,プリアンブルで以下を定義しておく。

\usepackage{listings,jlisting}
\lstset{%
  basicstyle={\ttfamily\normalsize},%
  frame={tbrl},
  numbers=left,%
  xleftmargin=3zw,%
  lineskip=-0.5ex%
}

これだけではだめで,

/usr/local/lib/python3.6/site-packages/pweave/formatters/tex.py

の,initformat の赤字部分を以下のように変更する。

    def initformat(self):
        self.formatdict = dict(codestart='\\begin{lstlisting}',
                               codeend='\\end{lstlisting}',
                               outputstart='\\begin{quote}\n\\begin{verbatim}',
                               outputend='\\end{verbatim}\n\\end{quote}',
                               figfmt='.pdf',
                               width='\\linewidth',
                               doctype='tex')

これで

<<>>=
x = 2
y = 4
print(x-y)
@

が,

のようなイメージで表示される。

すなわち,プログラム〔入力〕が行番号付きで,出力は行番号なしで出力される。

\color{} を使って,入力と出力を分けることもできるけど,色を使うのは避けたい。

コメント

Pweave のコードチャンク・オプション

2018年09月04日 | RとLaTeX
  • ネームなしのオプションは name または label
    すなわち,以下の 3 つ,name="aaa", label="aaa", aaa は同じもので,チャンクの名前を定義する
  • echo=True or (False)
    コードチャンクのエコー出力をするかどうか
  • evaluate=True or (False)
    コードチャンクを評価するかどうか
  • results="verbatim" or ("hidden")
    結果をテキスト出力("verbatim"),出力しない("hidden")
    "tex", "rst" は未実装?
  • term=False or (True)
    False のときは,結果はまとめて出力される
    True のときは,結果はその都度出力される
  • include=True or (False)
    True のときは,生成される図は自動的に文書に挿入される
    False のときは,図は生成されるが文書には挿入されない
  • fig=True or (False)
    True のときは,matplotlib plot で生成される図をファイルに書き出す
    このファイルは \includegraphics で文書に挿入できる
    figure 環境の caption は次の caption オプションで定義する
  • caption=""
    fig オプションと一緒に使い,figure 環境の caption を定義する
  • width
    文書に挿入される図の大きさ(任意の単位: cm, pt, linewidth など)
  • f_size=(8,6)
    保存される matplotlib の図の大きさ(単位はインチ)
    横幅,縦軸をタプルで指定
  • f_spines=True
    False のときは,matplot の図の右,上の spine を取り除く
  • f_env
    図の環境を追加する
  • f_pos="htbp"
    figure 環境の図の位置を指定する
  • wrap=True or (False, "code", "results")
    True ならば,長い行は 75 文字で折り返す
    "code","results" は指定されたものだけを折り返す
  • complete=True or (False)
  • source=ファイル名または Python モジュール・ファイル
    入力ソースをファイルやモジュール・ファイルにする

コメント

pythontex というのがあった

2018年09月01日 | ブログラミング

TeXShop の engines フォルダを見たら,pythontex というのがあった。

ちょっと,調べてみる。

なんとなれば,Psweave は 入力に対する出力はまとめられてしまうからなあ。

そのほかの不都合は前述の通り。

 

コメント

Pweave の TeXShop engine

2018年09月01日 | RとLaTeX

Sweave engine の最初の部分だけを以下のように変更して使う。

mp = function(file, makeindex=FALSE, silent=FALSE, deletePdfs=FALSE, deleteWorkfiles=FALSE, ...) {
  Sys.setlocale("LC_ALL", "ja_JP.UTF-8")
  if (grepl("\\.", file) == FALSE) {
    file = paste(file, "Rnw", sep=".")
  }
  cat("Input file:", file, "\n")
  cmd = sprintf("pweave -f tex %s", file)
  system(cmd)
  base = sub(".Rnw", "", file)
  :

 

コメント

Pweave --- Sweave for python

2018年09月01日 | RとLaTeX

R のための Sweave と同じような感じで,python のための Pweave(http://mpastell.com/pweave/) を使ってみる。

以下で,全角の <,> は,実際は半角で

TeXShop を使って,*.Pnw --> *.tex ---> *.pdf の操作を行うにあたって,大問題がある。

それは,TeXShop は .Pnw などという拡張子を認識しないということである(.Rnw は認識するのに)。

仕方ないから,.Rnw を流用することにするが,何とかならないか?

● チャンクの記述法

もう一つの問題点。

コマンドラインでは ```pythpn 〜 ``` が有効で,<<>>= 〜 @ はだめ。

TeXShop では全く逆。```pythpn 〜 ``` がだめで,<<>>= 〜 @ が有効。


次の 3 つの例は,だめ。

(1)

```python
a = 12
b = 34
print(a+b)
```

(2)

オプションを定義することもできる

```python, echo=False
a = 12
b = 34
print(a+b)
```

図を挿入する方法

```python
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 2*np.pi)
plt.plot(x, np.sin(x))
```

\bigskip

次の 2 つの例は有効。

<<fig = True, width = '12 cm', echo = True>>=
from pylab import *
plot(arange(10))
show()
@

<<>>=
x = 2
y = 4
print(x-y)
@


● インラインコード

コードをその場で評価する sqrt(8)= <% np.sqrt(8) %>
            8 の平方根を,この後へ挿入 sqrt(8)=2.8284271247461903


式を評価した結果を出力する sqrt(8)=<%= np.sqrt(8) %>
            8 の平方根を,この後へ挿入 sqrt(8)=2.8284271247461903

違いはないようにみえるのだが?

コメント

Pweave --- Sweave for python

2018年09月01日 | ブログラミング

R のための Sweave と同じような感じで,python のための Pweave(http://mpastell.com/pweave/) を使ってみる。

以下で,全角の <,> は,実際は半角で

TeXShop を使って,*.Pnw --> *.tex ---> *.pdf の操作を行うにあたって,大問題がある。

それは,TeXShop は .Pnw などという拡張子を認識しないということである(.Rnw は認識するのに)。

仕方ないから,.Rnw を流用することにするが,何とかならないか?

● チャンクの記述法

もう一つの問題点。

コマンドラインでは ```pythpn 〜 ``` が有効で,<<>>= 〜 @ はだめ。

TeXShop では全く逆。```pythpn 〜 ``` がだめで,<<>>= 〜 @ が有効。


次の 3 つの例は,だめ。

(1)

```python
a = 12
b = 34
print(a+b)
```

(2)

オプションを定義することもできる

```python, echo=False
a = 12
b = 34
print(a+b)
```

図を挿入する方法

```python
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 2*np.pi)
plt.plot(x, np.sin(x))
```

\bigskip

次の 2 つの例は有効。

<<fig = True, width = '12 cm', echo = True>>=
from pylab import *
plot(arange(10))
show()
@

<<>>=
x = 2
y = 4
print(x-y)
@


● インラインコード

コードをその場で評価する sqrt(8)= <% np.sqrt(8) %>
            8 の平方根を,この後へ挿入 sqrt(8)=2.8284271247461903


式を評価した結果を出力する sqrt(8)=<%= np.sqrt(8) %>
            8 の平方根を,この後へ挿入 sqrt(8)=2.8284271247461903

違いはないようにみえるのだが?

コメント

TeXShop で Sweave

2018年08月20日 | RとLaTeX

TeXShop で用意されている engine に Sweave 用のものがないので,作ってみた。

基本的には,以下のようであるが,肝は Sweave 関数の前に,Sys.setlocale("LC_ALL", "ja_JP.UTF-8") を置くこと。これがないと locale が "C" になっているので,日本語が正しく処理されない。

このファイルを Sweave.engine と名付けて,Engines フォルダに放り込んでおく。選択肢に "Sweave" が加わる。

#!/usr/local/bin/Rscript --vanilla --default-packages=methods,datasets,utils,grDevices,graphics,stats
file = sub("\\.Rnw", "", commandArgs(TRUE)[1])
Sys.setlocale("LC_ALL", "ja_JP.UTF-8")
Sweave(file, encoding="utf-8")
system(sprintf("/Library/TeX/texbin/uplatex %s.tex", file))
system(sprintf("/Library/TeX/texbin/dvipdfmx %s.dvi", file))

 

コメント

Sweave.sty

2018年07月29日 | RとLaTeX

Sweave で,

<<>>==
ans <- lm(Petal.Length ~ ., data=iris)
logLik(ans)
@

のようなチャンクを含む *.Rnw を

Sweave ==> uplatex ==> dvipdfmx すると,*.dvi ファイルを *.pdf ファイルにする際に

kpathsea: Running mktexpk --mfmode / --bdpi 600 --mag 1+0/600 --dpi 600 tcrm1000
mktexpk: Cannot find mktex.opt; check your installation.
kpathsea: Appending font creation commands to missfont.log.

dvipdfmx:warning: Could not locate a virtual/physical font for TFM "tcrm1000".
dvipdfmx:warning: >> There are no valid font mapping entry for this font.
dvipdfmx:warning: >> Font file name "tcrm1000" was assumed but failed to locate that font.
dvipdfmx:fatal: Cannot proceed without .vf or "physical" font for PDF output...

Output file removed.

なるエラーメッセージを吐いてとまってしまう。

R コンソールでは,直前のエラーメッセージが赤字で目立つのだが,その前に本当のエラーが暗示されている。

LaTeX Font Warning: Font shape `TS1/aett/m/n' undefined
(Font)              using `TS1/cmr/m/n' instead
(Font)              for symbol `textquotesingle' on input line 31.

ドツボにはまっていたのだが,`textquotesingle' だの,aett だのと書かれているのがポイント。

確かに,件のチャンクの出力は,

'log Lik.' -9.283269 (df=7)

と,シングルクオートを含む。さらに,そういえば Sweave.sty で ae が悪さをすることがあるというのを思い出した。

Sweave.sty の 8 行目,

\setboolean{Sweave@ae}{true}

\setboolean{Sweave@ae}{false}

にすることで,問題解決...したのかな?



コメント

文字列としての連結と"+"演算子

2018年06月24日 | ブログラミング

"+" を文字列連結演算子とする場合

"+" = function(e1, e2) {
  if (is.numeric(e1) && is.numeric(e2)) {
    base::"+"(e1, e2)
  } else {
    paste0(e1, e2)
  }
}

> "abc" + "12345" + "あいうえお"
[1] "abc12345あいうえお"

> 123456 + "numeric"
[1] "123456numeric"

> 456 + 100
[1] 556

> 456L + 100L
[1] 556

> 456 + 0
[1] 456

> TRUE + FALSE
[1] "TRUEFALSE"

> "123" + TRUE + "asd"
[1] "123TRUEasd"

> (1 == 1) + (2 != 3)
[1] "TRUETRUE"

 

コメント

文字列としての連結と"&"演算子

2018年06月10日 | ブログラミング

文字列と文字列の連結は,言語によって異なる記号が演算子として使われているが,R では二項演算子ではなく paste 関数,paste0 関数が使われる。

> paste("asd", "poi", sep="")
[1] "asdpoi"

> paste0("asd", "poi")
[1] "asdpoi"

> paste0(123, "asd")
[1] "123asd"

> paste0(123, TRUE)
[1] "123TRUE"

> paste0(123, T)
[1] "123TRUE"

> paste0(TRUE, FALSE)
[1] "TRUEFALSE"

TRUE と FALSE を連結するなんてことは普通は考えないこと。

& 記号を文字列としての連結演算子と定義する。ただし,2 項ともに論理型の場合には本来の & 演算子として使うというような関数を書いてみる。

"&" を "|" に変えても同じように動く。お好きな方を。

"&" = function(e1, e2) {
  if (is.logical(e1) && is.logical(e2)) {
    base::"&"(e1, e2)
  } else {
    paste0(e1, e2)
  }
}

> "abc" & "12345" & "あいうえお"
[1] "abc12345あいうえお"

> 123456 & "numeric"
[1] "123456numeric"

> 456 & 100
[1] "456100"

> 456 & 0
[1] "4560"

> TRUE & FALSE
[1] FALSE

> "123" & TRUE & "asd"
[1] "123TRUEasd"

> (1 == 1) & (2 != 3)
[1] TRUE

 

コメント

大きい数の積や商(その2)

2018年06月08日 | ブログラミング

超幾何分布


要因 a, b が独立な場合,以下のような分割表が得られる確率を求める

          B     not(B)   sum
A         x       o       m
not(A)  (k-x)     o       n
sum       k       o     (m+n)

超幾何分布により,求める確率は xCm * xCk-x / m+nCk

R の dhyper の引数はちょっと変で(?),x が生じる確率は dhyper(x, m, n, k)

          B     not(B)   sum
A         2       o       10
not(A)   (1)      o       5
sum       3       o      (15)

> x = 2; m = 10; n = 5; k = 3
> dhyper(x, m, n, k)
[1] 0.494505494505494
> choose(m, x) * choose(n, k-x) / choose(m+n, k)
[1] 0.494505494505495
> exp(lchoose(m, x) + lchoose(n, k-x) - lchoose(m+n, k))
[1] 0.494505494505495

しかし,

> x = 200; m = 2000; n = 5000; k = 600
> dhyper(x, m, n, k)
[1] 0.00104601828356522
> choose(m, x) * choose(n, k-x) / choose(m+n, k)
[1] NaN
> exp(lchoose(m, x) + lchoose(n, k-x) - lchoose(m+n, k))
[1] 0.00104601828356506

となり,普通のやり方では精度が足りない。

Rmpfr パッケージを用いて,単純に計算すると以下のようになる。

library(Rmpfr)
precBits=13
dhyper2 = function(x, m, n, k) {
  choose2 = function(a, b) {
    prod(mpfr(1:a, precBits)) / prod(mpfr(1:b, precBits)) / prod(mpfr(1:(a-b), precBits))
  }
  as.numeric(choose2(m, x) * choose2(n, k-x) / choose2(m+n, k))
}

> x = 200; m = 2000; n = 5000; k = 600
> dhyper(x, m, n, k)
[1] 0.00104601828356522
> dhyper2(x, m, n, k)
[1] 0.00104601828356522

dhyper2 の実行時間は 0.15 秒



コメント

大きい数の積や商

2018年06月07日 | ブログラミング

二項分布 B(x, n, p) = nCx p^x (1-p)^(n-x) の計算
R では,dbinom(x, n, p) で速く,簡単に求まる
これを,自分で計算してみるとどうなるかをやってみる。

式をそのまま記述すると以下のようになる。

options(digits = 15)

f = function(x, n, p) {
  return(choose(n, x) * p ^ x * (1 - p) ^ (n - x))
}

n が小さいうちは,これでも十分である。

> x = 0
> n = 10
> p = 0.5
> dbinom(x, n, p)
[1] 0.0009765625
> f(x, n, p)
[1] 0.0009765625

しかし,n が大きくなると,これじゃダメだ。

> x = 4000
> n = 10000
> p = 0.4
> dbinom(x, n, p)
[1] 0.00814316030659453
> f(x, n, p)
[1] NaN

==============================================

一つのやりかたは,掛算・割り算を log の足し算・引き算で計算し,結果の exp をとる。

g = function(x, n, p) {
  g0 = function(m) {
    ifelse(m <= 1, 0, sum(log(2:m))) # R の仕様からくる問題を避けるために関数を定義
  }
  return(exp(g0(n) - g0(x) - g0(n - x) + x * log(p) + (n - x) * log(1 - p)))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> g(x, n, p)
[1] 0.00814316030660582

このやり方では,有効数字は 10 桁ほどである。

==============================================

もう一つのやり方は,n の階乗は n+1 に対するΓ関数なので,それを利用して対数Γ関数を使う。

h = function(x, n, p) {
  return(exp(lgamma(n + 1) - lgamma(x + 1) - lgamma(n - x + 1) +
           x * log(p) + (n - x) * log(1 - p)))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> h(x, n, p)
[1] 0.00814316030654657

このやり方では,有効数字は 11 桁ほどである。

==============================================

R では nCx の対数は lchoose 関数で計算できるので,少し簡単に書ける。

v = function(x, n, p) {
  return(exp(lchoose(n, x) + x * log(p) + (n - x) * log(1 - p)))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> v(x, n, p)
[1] 0.00814316030660582

このやり方では,有効数字は 10 桁ほどである。

==============================================

nCx は整数値(約分すると分母は 1 になる)ということで,約分後の分子の log の和を求めて lchoose と同じ答えを目指す。

gcd = function(m, n) {
  while ((temp <- n %% m) != 0) {
    n <- m
    m <- temp
  }
  return(m)
}

lchoose2 = function(n, k) {
  k = min(k, n - k)
  if (k == 0) {
    return(0)
  }
  numerator = (n - k + 1):n
  denominator = 1:k
  for (den in seq_along(denominator)) {
    for (num in seq_along(numerator)) {
      GCD = gcd(denominator[den], numerator[num])
      if (GCD != 1) {
        denominator[den] = denominator[den] / GCD
        numerator[num] = numerator[num] / GCD
        if (denominator[den] == 1) {
          break
        }
      }
    }
    numerator = numerator[numerator != 1]
  }
  return(sum(log(numerator)))
}

w = function(x, n, p) {
  return(exp(lchoose2(n, x) + x * log(p) + (n - x) * log(1 - p)))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> w(x, n, p)
[1] 0.00814316030659842

このやり方では,有効数字は 12 桁に増加するが,計算時間が 5 秒近くかかる。
dbinom の結果の違いは,計算途中での計算精度が不足しているからである。

==============================================

R には多倍長演算をサポートするパッケージがいくつかあるが,そのうちの Rmfpr を使ってみる。
計算精度を適当(適切!)に指定して以下のような結果になる。

積・商を対数の和・差にする場合。

w2 = function(x, n, p) {
  library(Rmpfr)
  precBits = 70
  P = mpfr(p, precBits)
  lchoose = sum(log(mpfr((n-x+1):n, precBits))) - sum(log(mpfr(1:x, precBits)))
  as.numeric(exp(lchoose + x * log(P) + (n - x) * log(1 - P)))
}

積・商をそのまま計算する場合。

w3 = function(x, n, p) {
  library(Rmpfr)
  precBits = 70
  P = mpfr(p, precBits)
  choose = prod(mpfr((n-x+1):n, precBits)) / prod(mpfr(1:x, precBits))
  as.numeric(choose * (P ^ x) * (1 - P) ^ (n - x))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> w2(x, n, p)
[1] 0.00814316030659453
> w3(x, n, p)
[1] 0.00814316030659453

いずれでも, bdinom と同じ精度の結果が得られる。
計算時間は,積・商をそのまま計算する w3 の方が速く 0.14 秒ほどである(dbinom は瞬時)。

==============================================

二項分布 B(x, n, p) = nCx p^x (1-p)^(n-x) の計算の場合は dbinom があるので,以上のような事をやる必要はさらさらないが,「ものすごく大きな数とものすごく小さな数の積がほどほどの値になる場合」や「ものすごく大きな 2 つの数の商がほどほどの値になる場合」には,「積と商を対数の和と対数の差をとって結果の逆対数をとる」というのが必要になることもある。
そういうときには Rmfpr を使うのがもっとも効果的のようである。
... と,そういうこと。

コメント