modulus は「絶対値」じゃない

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

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

> 30284005485024840 %% 7
[1] 4

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

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

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

いくら，例だからと言って，こんな汚い例はない方がよい。

> 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))

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

> 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;
}

百五減算

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(Thu) 17:48 No. 13352

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

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

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

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

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

平均値だけが同じだの・違うだのいっても，しようがないと思いますけど？

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

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 曲線

前に作ったはずなのに見つからなかった
ドラゴンカーブのプログラムとほんのちょっとの違いで，結果はかなり違うものになる

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

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

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

フラクタルの例として木

これも以前書いたもの

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)
}

描画結果

コッホ曲線

これも，以前書いたもの

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))

描画結果

ドラゴン曲線

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)

描画結果

スピログラフ（ついでに）

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)

描画結果

ユークリッドの互除法

良くあるパターンであるが，引数を入れ替える必要はない

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 を使うべき場合

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