裏 RjpWiki

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

modulus は「絶対値」じゃない

2010年09月09日 | ブログラミング
絶対値で精度が完全に失われた可能性があります 」というのは実に日本語的に意味が分からない珍訳だと思うのだけど,英語では "probable complete loss of accuracy in modulus" だ。
modulus が絶対値という訳になるのは,「(複素数の)絶対値」だから,当てはまらないだろう。
modulus ということで,一番最初に頭に浮かぶのは「(整数論の)法」

それにしても,どういう状態か分からないのでいろいろ捜索の結果,以下のような場合に出たり出なかったりというエラーのようだ。

> 30284005485024840 %% 7
[1] 4

> 30284005485024840 %% 6
[1] 0
 警告メッセージ:
 絶対値で精度が完全に失われた可能性があります 

やはり,modulo に関係するもののようだけど,除数により出たり出なかったりというのがわからない。
コメント

処理に適したデータ構造を

2010年09月08日 | 裏 RjpWiki
いくら,例だからと言って,こんな汚い例はない方がよい。

> hoge1 <- 1; hoge2 <- 2; hoge3 <- 3; hoge4 <- 4
> N <- 4
> X <- vector("list",N)
> for (i in 1:N) X[[i]] <- eval(parse(text=paste("hoge",i,sep="")))
> str(X)
List of 4
 $ : num 1
 $ : num 2
 $ : num 3
 $ : num 4
> foo <- function(x) x^2
> Y <- lapply(X,foo)
> str(Y)
List of 4
 $ : num 1
 $ : num 4
 $ : num 9
 $ : num 16
> Y[[1]];Y[[2]];Y[[3]];Y[[4]]
[1] 1
[1] 4
[1] 9
[1] 16


X はベクトルでよいし,結果も lapply でなく sapply でベクトルに付値するべし

hoge1 <- 1; hoge2 <- 2; hoge3 <- 3; hoge4 <- 4
N <- 4
for (i in 1:N) X[i] <- eval(parse(text=paste("hoge",i,sep="")))
foo <- function(x) x^2
(Y <- sapply(X,foo))
コメント

計算誤差とはそういうもの

2010年09月03日 | ブログラミング
> options(digits=20)
> (x <- sum(1/factorial(0:18)))
[1] 2.718281828459045
> x == exp(1)
[1] TRUE
> (y <- 1+sum(1/factorial(1:18)))
[1] 2.718281828459046
> y == exp(1)
[1] FALSE

ちなみに,計算精度の範囲内では,

> (x <- sum(1/factorial(0:17)))
[1] 2.718281828459045
> x == exp(1)
[1] TRUE

で,十分

ちなみに,奥村先生のアルゴリズム入門にある C プログラムは長い。
数学の教科書に書かれていることを確認してみようかなあと思う人も,C で確かめようとすると大変だけど,R だと超簡単。

/***********************************************************
    e.c -- 自然対数の底
***********************************************************/

long double ee(void)
{
    int n;
    long double e, a, prev;

    e = 0;  a = 1;  n = 1;
    do {
        prev = e;  e += a;  a /= n;  n++;
    } while (e != prev);
    return e;
}

#include <stdio.h>
#include <stdlib.h>
#include <float.h>

int main()
{
    printf("e = %.*Lg¥n", LDBL_DIG, ee());
    return EXIT_SUCCESS;
}
コメント

百五減算

2010年09月02日 | ブログラミング
1 から 100 までの整数を考えよう。
その数を,3,5,7 で割った余り a, b, c を
百五減算(a, b, c) と入力すると,考えた整数が判明する。

> # 関数の定義
> 百五減算 <- function(a, b, c) (70*a+21*b+15*c) %% 105

> 百五減算(2, 3, 2)
[1] 23

> 百五減算(1, 4, 0)
[1] 49
コメント

分布を識別するのは,平均値だけじゃないからなあ

2010年09月02日 | 統計学
1以下の数の期待値に対する検定方法   
投稿者:藤田 2010/09/02(Thu) 17:48 No. 13352

    0から1の間の値(X軸)をとるべき分布形状(最大値があるのでべき分布といえるか不明のため)の確率分布A(分布自体は不明)があり。
    この期待値μ1が分っているとします。
    ここでサンプル値が100個(x1~x100)得られたとき、これらがAから得られたことを帰無仮説とする検定をしたいです。
    期待値があるのでカイ2乗検定を考えましたが、小数となるためうまく使えそうにありません。

    恐れ入りますが、お知恵をいただけると幸いです

Re: 1以下の数の期待値に対する検定方法
投稿者:ひの 2010/09/02(Thu) 18:45 No. 13353

    母平均の検定(母分散未知の場合)ですね。

    http://aoki2.si.gunma-u.ac.jp/lecture/Average/Mean1.html


平均値だけが同じだの・違うだのいっても,しようがないと思いますけど?
コメント

複数の関数は使う順序も考えよう

2010年09月02日 | ブログラミング
char <- c(rep("a", 2), rep("b", 3), rep("c", 4))    # 関数は入れ子にすることができる!
char
#=> [1] "a" "a" "b" "b" "b" "c" "c" "c" "c"

以下に示すような方法を使うのがよい

> rep(c("a", "b", "c"), c(2, 3, 4))
[1] "a" "a" "b" "b" "b" "c" "c" "c" "c"

この場合には letters を使うこともできるけど
> rep(letters[1:3], c(2, 3, 4))
[1] "a" "a" "b" "b" "b" "c" "c" "c" "c"

例に挙げられたくらいなら良いけど,
rep(letters, 1:26)
とか
rep(letters, each=3)
のようなものを例に挙げたような方法で書こうとすると,音を上げるだろう
コメント

C 曲線

2010年09月01日 | ブログラミング
前に作ったはずなのに見つからなかった
ドラゴンカーブのプログラムとほんのちょっとの違いで,結果はかなり違うものになる

drawC <- function(a.x, a.y, b.x, b.y, n)
{
    x <- b.x-a.x
    y <- a.y-b.y
    c.x <- a.x+(x-y)/2
    c.y <- b.y-(x-y)/2
    if (n == 0) {
        lines(c(a.x, c.x, b.x), c(a.y, c.y, b.y), col="red")
    }
    else {
        drawC(a.x, a.y, c.x, c.y, n-1)
        drawC(c.x, c.y, b.x, b.y, n-1)
    }
}
old <- par(mar=c(0, 0, 0, 0))
plot(c(-6, 4), c(-8, 8), type="n", axes=FALSE, xlab="", ylab="", asp=1)
drawC(0, 4, 0, -4, 15)
par(old)

描画結果

n=15

n=12

n=10
コメント

πを求めるシミュレーション

2010年09月01日 | ブログラミング
Taglibro de H
モンテカルロ法
http://ito-hi.blog.so-net.ne.jp/2006-09-05


短くすればよいと言うものではないが,二次元配列を利用すると良い
配列の二乗,列和に colSums,列和が 1 以下になる平均値(!)を求め(sum(p)/n って,mean(p) のこと)4倍する

set.seed(124)
n <- 100000
mean(colSums(matrix(runif(2*n), 2)^2) < 1)*4

答え 3.14484
コメント

フラクタルの例として木

2010年09月01日 | ブログラミング
これも以前書いたもの

drawTree <- function(a.x, a.y, b.x, b.y, n, col=1)
{
    STEM.RATIO <- 0.25
    BRANCH.RATIO <- 0.6
    xx <- b.x-a.x
    yy <- a.y-b.y
    angle1 <- atan(yy/xx)+pi/6
    angle2 <- atan(yy/xx)-pi/6*1.5
    center.length <- sqrt(xx^2+yy^2)*(1-STEM.RATIO)
    branch.length <- BRANCH.RATIO*center.length
    sig <- ifelse(xx >= 0, 1, -1)
    c.x <- a.x+STEM.RATIO*xx
    c.y <- a.y-STEM.RATIO*yy
    d.x <- c.x+sig*(branch.length*cos(angle1))
    d.y <- c.y-sig*(branch.length*sin(angle1))
    e.x <- c.x+sig*(branch.length*cos(angle2))
    e.y <- c.y-sig*(branch.length*sin(angle2))
    segments(a.x, a.y, c.x, c.y, col=col)
    if (n <= 0) {
        segments(c.x, c.y, c(b.x, d.x, e.x), c(b.y, d.y, e.y), col=col)
    }
    else {
        drawTree(c.x, c.y, b.x, b.y, n-1, col=col)
        drawTree(c.x, c.y, d.x, d.y, n-1, col=col)
        drawTree(c.x, c.y, e.x, e.y, n-1, col=col)
    }
}
plot(c(00,500), c(0, 500), type="n", axes=FALSE, xlab="", ylab="", asp=1)
for (i in 1:10) {
    x <- sample(500, 1)
    y1 <- sample(500, 1)
    y2 <- y1-runif(1, min=100, max=400)
    drawTree(x, y2, x, y1, sample(4, 1)+2, col=i)
}

描画結果

コメント

コッホ曲線

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

これも,以前書いたもの

drawKoch <- function(a.x, a.y, b.x, b.y, n)
{
    c.x <- (2*a.x+b.x)/3
    c.y <- (2*a.y+b.y)/3
    d.x <- (a.x+2*b.x)/3
    d.y <- (a.y+2*b.y)/3
    x <- b.x-a.x
    y <- a.y-b.y
    d <- sqrt(x^2+y^2)/sqrt(3)
    if (x >= 0) {
        a1 <- atan(y/x)+pi/6
        e.x <- a.x+d*cos(a1)
        e.y <- a.y-d*sin(a1)
    }
    else {
        a2 <- atan(y/x)-pi/6
        e.x <- b.x+d*cos(a2)
        e.y <- b.y-d*sin(a2)
    }
    if (n <= 0) {
        lines(c(a.x, c.x, e.x, d.x, b.x), c(a.y, c.y, e.y, d.y, b.y), col="red")
    }
    else {
        drawKoch(a.x, a.y, c.x, c.y, n-1)
        drawKoch(c.x, c.y, e.x, e.y, n-1)
        drawKoch(e.x, e.y, d.x, d.y, n-1)
        drawKoch(d.x, d.y, b.x, b.y, n-1)
    }
}
plot(c(50, 450), c(70, 428), type="n", axes=FALSE, xlab="", ylab="", asp=1)
r <- 180
ox <- oy <- 250
p.x <- r*cos(pi*(3/6))+ox
p.y <- r*sin(pi*(3/6))+oy
q.x <- r*cos(pi*(7/6))+ox
q.y <- r*sin(pi*(7/6))+oy
r.x <- r*cos(pi*(11/6))+ox
r.y <- r*sin(pi*(11/6))+oy
n <- 7
drawKoch(p.x, p.y, q.x, q.y, n)
drawKoch(q.x, q.y, r.x, r.y, n)
drawKoch(r.x, r.y, p.x, p.y, n)
grid(40)
points(c(p.x, q.x, r.x), c(p.y, q.y, r.y))

描画結果

コメント

ドラゴン曲線

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

Taglibro de H
ドラゴン曲線
http://ito-hi.blog.so-net.ne.jp/2008-09-06


前に作っていたものだけど

drawDragon <- function(a.x, a.y, b.x, b.y, n)
{
    x <- b.x-a.x
    y <- a.y-b.y
    c.x <- a.x+(x+y)/2
    c.y <- b.y+(x+y)/2
    if (n <= 0) {
        lines(c(a.x, c.x, b.x), c(a.y, c.y, b.y), col="red")
    }
    else {
        drawDragon(a.x, a.y, c.x, c.y, n-1)
        drawDragon(b.x, b.y, c.x, c.y, n-1)
    }
}
plot(c(50, 450), c(70, 428), type="n", axes=FALSE, xlab="", ylab="", asp=1)
p.x <- 170
p.y <- 140
q.x <- 400
q.y <- 350
drawDragon(p.x, p.y, q.x, q.y, 10)
grid(40)

描画結果

コメント

スピログラフ(ついでに)

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

Taglibro de H
http://ito-hi.blog.so-net.ne.jp/2009-11-05

ついでにスピログラフのプログラムを while などを使わずに書いてみる。
ずいぶん短くなるし,描画もあっという間であるけれど。
while で描く方が,とろとろ描いていておもしろいのだけど(笑)。

spirograph <- function(
    gr,                # 外円の半径
    sr,                # 内円の半径
    pr,                # 内円内のペンの位置
    accuracy=100,        # 描画精度
    col = "black",        # 描画色
    lwd = 1)            # 描画線幅
{
    gcd <- function(a, b) {
        r <- a %% b
        return(ifelse(r == 0, b, gcd(b, r)))
    }

    gr <- trunc(gr)
    sr <- trunc(sr)
    pr <- trunc(pr)
    if (gr <= 0 | sr <= 0 | pr <= 0) {
        stop("Parameters must be postive.")
    }
    theta0 <- seq(0, 2*pi*sr/gcd(gr, sr), length=accuracy*sr/gcd(gr, sr))
    theta1 <- theta0*(gr/sr)
    x <- (gr-sr)*sin(theta0)-pr*sin(theta1)
    y <- (sr-gr)*cos(theta0)-pr*cos(theta1)
    old <- par(mar=c(0, 0, 0, 0))
    plot(x, y, type="n", axes=FALSE, xlab="", ylab="", asp=1)
    lines(x, y, lwd=lwd, col=col)
    par(old)
}

layout(matrix(1:4, 2))
spirograph(150, 108, 60)
spirograph(72, 552, 600, lwd = 3)
spirograph(11, 30, 45, col="red")
spirograph(15*29, 15*23, 50)
layout(1)

描画結果

コメント

ユークリッドの互除法

2010年09月01日 | ブログラミング
良くあるパターンであるが,引数を入れ替える必要はない

Taglibro de H
http://ito-hi.blog.so-net.ne.jp/2009-11-05

  ## 最大公約数
  gcd <- function(a, b) {
    if (a < b) {
      t <- a
      a <- b
      b <- t
    }
    r <- a %% b
    return(ifelse(r == 0, b, gcd(b, r)))
  }


a が b より小さい場合でも,一回余分に関数呼び出しがあるだけで,ちゃんと解が求まる。

  gcd2 <- function(a, b) { # これでいいのだ
    r <- a %% b
    return(ifelse(r == 0, b, gcd2(b, r)))
  }


> gcd(168, 48)
[1] 24
> gcd(48, 168)
[1] 24

> gcd2(168, 48)
[1] 24
> gcd2(48, 168)
[1] 24
コメント

matrix ではなく data.frame を使うべき場合

2010年09月01日 | ブログラミング
Taglibro de H

http://ito-hi.blog.so-net.ne.jp/2008-12-02?comment_success=2010-09-01T15:26:42&time=1283322402
R: order()で並べ替え [統計] [編集]

こういうデータがあったとする。

m <- matrix(c(
    261, "コナラ","Quercus serrata",
     65, "ウリハダカエデ","Acer rufinerve",
    288, "アケビ","Akebia quinata",
    177, "モチツツジ","Rhododendron macrosepalum",
    288, "ミツバアケビ","Akebia trifoliata",
     51, "タラノキ","Aralia elata",
    261, "アラカシ","Quercus glauca",
    261, "アベマキ","Quercus variabilis",
    177, "コバノミツバツツジ","Rhododendron reticulatum",
    171, "エゴノキ","Styrax japonica"),
  ncol = 3, byrow = TRUE)

途中省略

ただし、日本語ではうまく並び替えできない。
> m[order(as.numeric(m[,1]), m[,2]),]
      [,1]  [,2]                 [,3]                      
 [1,] "51"  "タラノキ"           "Aralia elata"            
 [2,] "65"  "ウリハダカエデ"     "Acer rufinerve"          
 [3,] "171" "エゴノキ"           "Styrax japonica"         
 [4,] "177" "モチツツジ"         "Rhododendron macrosepalum"
 [5,] "177" "コバノミツバツツジ" "Rhododendron reticulatum"
 [6,] "261" "コナラ"             "Quercus serrata"         
 [7,] "261" "アベマキ"           "Quercus variabilis"      
 [8,] "261" "アラカシ"           "Quercus glauca"          
 [9,] "288" "アケビ"             "Akebia quinata"          
[10,] "288" "ミツバアケビ"       "Akebia trifoliata"  


matrix では 1 列目も character になってしまっている。
また,2 列目は character だけど,エンコーディングのバイトが辞書順になっていない。
解決法は,元のデータを data.frame として定義する。そうすれば,うまくいく。

> m <- read.csv(stdin(), header=FALSE)
0:     261, "コナラ","Quercus serrata"
1:      65, "ウリハダカエデ","Acer rufinerve"
2:     288, "アケビ","Akebia quinata"
3:     177, "モチツツジ","Rhododendron macrosepalum"
4:     288, "ミツバアケビ","Akebia trifoliata"
6:      51, "タラノキ","Aralia elata"
7:     261, "アラカシ","Quercus glauca"
8:     261, "アベマキ","Quercus variabilis"
9:     177, "コバノミツバツツジ","Rhododendron reticulatum"
10:     171, "エゴノキ","Styrax japonica"
11:
>

> m
    V1                  V2                        V3
1  261              コナラ           Quercus serrata
2   65      ウリハダカエデ            Acer rufinerve
3  288              アケビ            Akebia quinata
4  177          モチツツジ Rhododendron macrosepalum
5  288        ミツバアケビ         Akebia trifoliata
6   51            タラノキ              Aralia elata
7  261            アラカシ            Quercus glauca
8  261            アベマキ        Quercus variabilis
9  177  コバノミツバツツジ  Rhododendron reticulatum
10 171            エゴノキ           Styrax japonica

> m[order(m[,1], m[,2]),]
    V1                  V2                        V3
6   51            タラノキ              Aralia elata
2   65      ウリハダカエデ            Acer rufinerve
10 171            エゴノキ           Styrax japonica
9  177  コバノミツバツツジ  Rhododendron reticulatum
4  177          モチツツジ Rhododendron macrosepalum
8  261            アベマキ        Quercus variabilis
7  261            アラカシ            Quercus glauca
1  261              コナラ           Quercus serrata
3  288              アケビ            Akebia quinata
5  288        ミツバアケビ         Akebia trifoliata

コメント