裏 RjpWiki

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

多倍長精度計算パッケージ Rmpfr

2014年07月30日 | ブログラミング

http://minato.sip21c.org/seek_pi_by_Archimedes.R
中澤さんの,「100ビット実数を使って円周率を10桁正しく出すには正何角形まで必要かを求めるためのRコード」について

Rmpfr というパッケージの使い方がこれだけでわかり,ありがたかった。

もともと,計算時間がかかるというプログラムではなかったが,ちょっとした工夫で計算時間が短くなることがわかったので,記録しておこう。

中澤さんのオリジナルプログラム

func1 = function() {
    library(Rmpfr)
    n = 6
    p = sqrt(mpfr(3, 100)) / mpfr(2, 100)
    a = sqrt(mpfr(1, 100) - p ^ 2)
    ad = a / p
    L = a * n
    M = ad * n
    while ((M - L) > mpfr(1e-10, 100)) {
        n = n * mpfr(2, 100)
        p = sqrt((mpfr(1, 100) + p) / mpfr(2, 100))
        a = sqrt(mpfr(1, 100) - p ^ mpfr(2, 100))
        ad = a / p
        L = a * n
        M = ad * n
    }
    invisible(c(n, L, M))
}

毎回 mpfr を使うのを避ける

func2 = function() {
    library(Rmpfr)
    n = mpfr(6, 100)
    one = mpfr(1, 100)
    two = mpfr(2, 100)
    epsilon = mpfr(1e-10, 100)
    p = sqrt(mpfr(3, 100)) / two
    a = sqrt(one - p ^ 2)
    ad = a / p
    L = a * n
    M = ad * n
    while ((M - L) > epsilon) {
        n = n * two
        p = sqrt((one + p) / two)
        a = sqrt(one - p ^ two)
        ad = a / p
        L = a * n
        M = ad * n
    }
    invisible(c(n, L, M))
}

混合演算の場合は,自動的に mpfr してくれるようなので,そのまま定数を書く

func3 = function() {
    library(Rmpfr)
    n = 6
    epsilon = mpfr(1e-10, 100)
    p = sqrt(mpfr(3, 100)) / 2
    a = sqrt(1 - p ^ 2)
    ad = a / p
    L = a * n
    M = ad * n
    while ((M - L) > epsilon) {
        n = n * 2
        p = sqrt((1 + p) / 2)
        a = sqrt(1 - p ^ 2)
        ad = a / p
        L = a * n
        M = ad * n
    }
    invisible(c(n, L, M))
}

更に,重箱の隅をつつき,計算時間を搾り取る

func4 = function() {
    library(Rmpfr)
    n = 6
    epsilon = mpfr(1e-10, 100)
    p = sqrt(mpfr(3, 100)) / 2
    a = sqrt(1 - p ^ 2)
    ad = a/p
    L = a * n
    M = ad * n
    while ((M - L) > epsilon) {
        n = n + n
        p = (1 + p) * 0.5
        a = sqrt(1 - p)
        p = sqrt(p)
        L = a * n
        M = L / p
    }
    invisible(c(n, L, M))
}

library(rbenchmark)
benchmark(func1(), func2(), func3(), func4())

結局の所,mpfr が不要な定数はそのまま書くのが速いらしい。

> benchmark(func1(), func2(), func3(), func4())
     test replications elapsed relative user.self sys.self
1 func1()          100  22.560    2.429    22.550    0.146
2 func2()          100  13.763    1.482    13.748    0.075
3 func3()          100  12.017    1.294    12.003    0.066
4 func4()          100   9.287    1.000     9.277    0.051

結果が正しいことを確認

> print(func1())
3 'mpfr' numbers of precision  100   bits
[1]                            786432 3.1415926535814376697324981306951 3.1415926536065043758325431868181
> print(func2())
3 'mpfr' numbers of precision  100   bits
[1]                            786432 3.1415926535814376697324981306951 3.1415926536065043758325431868181
> print(func3())
[[1]]
[1] 786432

[[2]]
'mpfr1' 3.1415926535814376697324981306951

[[3]]
'mpfr1' 3.1415926536065043758325431868181

> print(func4())
[[1]]
[1] 786432

[[2]]
'mpfr1' 3.1415926535814376697324981306951

[[3]]
'mpfr1' 3.1415926536065043758325431868181

コメント

Excel だって,賢く使おうョ

2014年07月29日 | ブログラミング

R で相関行列を(一気に)計算する
http://pepper.is.sci.toho-u.ac.jp/index.php?%A5%CE%A1%BC%A5%C8%2FR%2F%C1%EA%B4%D8%B9%D4%CE%F3cor
の前段で,Excel を使って相関係数行列を計算するやり方を,
7.4 相関行列を計算する(2)
http://kogolab.chillout.jp/elearn/icecream/chap7/sec4.html
を参照して書いているのだけど...

いくら Excel だからといって,PEARSON 関数を使って(関数式のコピーをするので,一つのセルごとに関数式を入れなくて済むよ!便利だね!!)なんて,ださい。

行列計算を使ってやった方が絶対便利。
テクニックは,データ領域に名前を付けること,関数の引数でその名前を使うこと,当然ながら,行列計算をやること。

あ,そうそう,分析ツールなんて使わないから(だって,Mac 用の Excel にはないから)

コメント

学術用語

2014年07月18日 | 統計学

> 統計モデリングのための堅牢な(または「耐性」)の方法は

という表現が某所であったが,robust のことなんだろうね。

通常は,「頑健な」という用語が使われると思う。

robustness は 頑健性

まあ,意味が通じればいいのではあるが。

コメント

named vector の活用

2014年07月14日 | ブログラミング

http://d.hatena.ne.jp/a_bicky/20140713/1405211663
R の named vector から subset を得る際の最悪計算量は O(n^2)

の前半部分で,データフレームの factor 列の情報に基づき新たな列を算出するのに,names を使うと速いよと書いてあった。

> nrow(students)
[1] 3000000
> class(students$grade)
[1] "factor"
> levels(students$grade)
[1] "可"   "不可" "優"   "良"  
> scores
  優   良   可 不可
  30   20   10    0

という状況で,students$grade が,優,良,可,不可 のとき,students$score が 30, 20, 10, 0 となるようにするということ。

> system.time(students$score <- scores[students$grade])
   ユーザ   システム       経過  
     0.087      0.002      0.089 


今まで,私なら,以下のようにしていた。

> system.time(students$score2 <- c(10, 0 , 30, 20)[as.integer(students$grade)])
   ユーザ   システム       経過  
     0.039      0.001      0.039


c(10, 0 , 30, 20) の順序は levels(students$grade) の順序に依存しているので,並び方がおかしいのはこのやりかたのせいではない。

倍速になったとはいっても,300 万行においてなので,意味のある差ではない。

scores[students$grade] の方がわかりやすいであろう。

コメント

素数の個数

2014年07月10日 | ブログラミング

某所で 2, 5, 10, 19, 54, 224, 312, 616, 888, 977 より小さい数字の中で,素数が幾つあるかを求める Java プログラムを募集している。
Java で書けといわれて書けるかどうかはプログラマの技量を測るには適当か(適切ではない)と思うが,こんな計算量の小さな問題を解くなら,R の gmp を使えば朝飯前。
何を使おうと,最短時間で解答を得られれば,何の問題も無いと思うけど。
で,Java での解を募集している訳なので,これをこのまま Java に出来るわけも無い(うん?Java からも gmp ライブラリを参照できるかな?)ので,ネタバレにはならないだろう。
「より小さい数字の中で」とあるので,ちょっとまぬけな対応をした。
> library(gmp)
> n = c(2, 5, 10, 19, 54, 224, 312, 616, 888, 977)
> sapply(n, function(m) sum(isprime(2:(m-1)) == 2))
 [1]   1   2   4   7  16  48  64 112 154 164
この答は,間違っている。2 より小さい素数は 0 個なのに 1 になっているから。
2 を特別扱いするのも,なんだかなあ...である
> sapply(n, function(m) if (m == 2) 0 else sum(isprime(2:(m-1)) == 2))
 [1]   0   2   4   7  16  48  64 112 154 164

 コメントをいただいた。最初に n-1 しなさいと。

そうですよね。気づかなかった。

> sapply(n-1, function(m) sum(isprime(1:m) == 2))
 [1]   0   2   4   7  16  48  64 112 154 164

コメント (1)

小ネタ:ベクトルから最後の k 個を削る

2014年07月09日 | ブログラミング

x = 1:20
k = 5

普通は,ベクトルの長さを求めてから必要部分を取り出す

n = length(x)
(y = x[1:(n-k)])

以下のようにすると,長さを求めなくてもよいし,記述もスマートか?

(y = head(x, -k))

tail を使うと以下のようなこともできる

(z = tail(x, -k))

コメント

CodeIQ 7 を数える。の解答例

2014年07月08日 | ブログラミング

CodeIQ 7 を数える
http://blog.goo.ne.jp/r-de-r/e/5de4a4c285dfdbf0dc86dedf85ebfcfd

ソースをどこに置いたか,失念していたのでちょっと遅れた。
基本的には,このようなプログラムは,ブルートフォースはもってのほか。
理論的に証明は難しくても,規則性を発見してそれを利用しなくては,有限の計算時間で解を求めるようなプログラムを得るのは難しい。

偉そうに言っていて,外していたら,非常に恥ずかしいかも知れないが,別に私は天才的なプログラマでもないし,ほかの人から「こんなことに気づかない程の馬鹿と思われても一向に気にしない(本当は,すごく気にして,落ち込んで,再起不能になる)。

で,解説抜きでプログラム提示(解説した文書も行方不明なんですよね...)

func1 = function(n) {
    return(sum(unlist(strsplit(as.character(1:n), "")) == "7"))
}

func2 = function(n) {
    ans = 0
    repeat {
        keta = nchar(n) - 1
        m = n %/% (10 ^ keta)
        ans = ans + m * keta * 10 ^ (keta - 1)
        if (m == 7) {
            ans = ans + n - 7 * 10 ^ keta + 1
        } else if (m > 7) {
            ans = ans + 10 ^ keta
        }
        n = n %% (10 ^ keta)
        if (n == 0)
            break
    }
    return(ans)
}

func1 は真っ正直なプログラムで,間違いない答えを出してくれる(と,思うし,思いたい)。けれども,計算時間がかかりすぎる。
func2 が,回答候補プログラム。

両方のプログラムで,実行可能範囲で両者の解が一致することを確認した。

それでもなお間違えていたとすれば,申し訳ない。

以上で~~す

コメント