裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

計算精度

2012年06月30日 | ブログラミング

一般的な電卓の計算精度というか扱える数値の最大値は我々が思っているよりは小さい。一般的には 1e100 程度だ。

コンピュータはそれよりは大きいものの限度はある。.Machine によれば,

> .Machine
  :
$double.xmax
[1] 1.797693e+308

ということなので,例えば,階乗は,170! までしか求められない。

> factorial(170)
[1] 7.257416e+306
> factorial(171)
[1] Inf
 警告メッセージ:
In factorial(171) :  'gammafn' 中の値が範囲を超えています

要するに,必要な計算をする場合に,途中でこの範囲(限界)を超えてはならないということだ。

例えば,2×2分割表の独立性を検定するときに,χ自乗統計量を計算 n*(ad-bc)/(a+b)/(c+d)/(a+c)/(b+d) するわけだが(記号法の説明は省略),その計算順序はどうであってもよいわけではない(普通はどうでもよいのだけど)。

> func1 <- function(x) {
+     y <- addmargins(x)
+     a  <- y[1,1]; b  <- y[1,2]; ab <- y[1,3]
+     c  <- y[2,1]; d  <- y[2,2]; cd <- y[2,3]
+     ac <- y[3,1]; bd <- y[3,2]; n  <- y[3,3]
+     return(n*(a*d-b*c)^2/(ab*cd*ac*bd))
+ }

> func2 <- function(x) {
+     y <- addmargins(x)
+     a  <- y[1,1]; b  <- y[1,2]; ab <- y[1,3]
+     c  <- y[2,1]; d  <- y[2,2]; cd <- y[2,3]
+     ac <- y[3,1]; bd <- y[3,2]; n  <- y[3,3]
+     return(n*(a*d-b*c)^2/ab/cd/ac/bd)
+ }

> func3 <- function(x) {
+     y <- addmargins(x)
+     a  <- y[1,1]; b  <- y[1,2]; ab <- y[1,3]
+     c  <- y[2,1]; d  <- y[2,2]; cd <- y[2,3]
+     ac <- y[3,1]; bd <- y[3,2]; n  <- y[3,3]
+     return(n/ab*(a*d-b*c)/cd/ac*(a*d-b*c)/bd)
+ }

普通の範囲ではどれでもちゃんと計算できる # 1-1, # 1-2, # 1-3

> x <- matrix(1:4, 2)
> func1(x) # 1-1
[1] 0.07936508
> func2(x) # 1-2
[1] 0.07936508
> func3(x) # 1-3
[1] 0.07936508

>しかし,例外的な場合には,計算順序を考えないと計算は失敗する # 2-1, # 2-2 は失敗する。しかし,# 2-3 はちゃんと計算できる。

x <- x*1e62
> func1(x) # 2-1
[1] Inf # 何てっこった
> func2(x) # 2-2
[1] Inf # 何てっこった
> func3(x) # 2-3
[1] 7.936508e+60 # 正しい

途中の計算を log を使って,一番最後に exp を使ってという方策で対処できる場合もあるけど,(アルゴリズムがマズイ場合は問題外だが),それでもダメ(途中の結果が限界を超える)ということもある。

「そのような場合は考慮外」とすましていることもできるけど,それは,やはり恥ずかしいことだよ。

どんな場合にも正しい答えを与えるようにプログラムする事は大変難しい(場合もある。プログラマのスキルにもよる)。

しかし,可能な限り,どんな場合にも正しい答えを出せるようにするのが「職人」だ。悔しければ,やってみな。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

メモリと計算時間のトレードオフ

2012年06月29日 | ブログラミング

DNA鑑定とバースデイ・パラドクス」にて

・対数を使うにしても、Rでやるにもn!;n=10^9とかになってくると、メモリの問題が出てくるようだ
・そんなのを少し回避しながら計算してみた

メモリと計算時間は互いにトレードオフの関係にある。
計算時間は掛かる(べらぼうに!でも許容範囲)がメモリは全く使わないプログラムを書けばよいだけだ。
このプログラムは引用元のプログラムが以下に無駄な計算をしているかを明示している。

func <- function(n) {
    q <- 4.7e12
    p <- 1
    for (i in 1:n - 1) {
        p <- p*(q-i)/q
    }

}

n = 10^i で i が 1 大きくなると,計算時間は 10 倍になる。
システム時間がほとんど掛かっていないところにも注目。

> for (i in 0:8) {
+     print(system.time(func(10^i)))
+ }
   ユーザ   システム       経過  
     0.001      0.000      0.000
   ユーザ   システム       経過  
         0          0          0
   ユーザ   システム       経過  
     0.001      0.000      0.001
   ユーザ   システム       経過  
     0.002      0.000      0.002
   ユーザ   システム       経過  
     0.018      0.000      0.018
   ユーザ   システム       経過  
     0.190      0.001      0.190
   ユーザ   システム       経過  
     1.941      0.011      1.939
   ユーザ   システム       経過  
    19.492      0.089     19.387
   ユーザ   システム       経過  
   196.289      1.623    196.516

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:ベクトル計算をしよう

2012年06月29日 | ブログラミング

DNA鑑定とバースデイ・パラドクス」にて

ns を 10^(seq(from=0,to=8,by=0.1)) に対して確率を求めているが,各 ns[i] について,毎回 1 からns[i] まで重複して計算をしている始めることになっている。無駄だ。
何を計算しているのか,よく考えよう。

> system.time({
+ # 4.7兆分の1の確率、と考えるための数字
+ q<-4.7*10^12
+
+ # 調べる人数(世界人口65億人だったり、日本人口1.2億人だったり、母集団100万人、10万人、1万人…など)
+ ns<- round(10^(seq(from=5,to=8,by=0.1)))
+
+ birthday<-rep(0,length(ns))
+ # 対数で計算しよう
+ # 同じ型のペアが少なくとも1組は存在する確率
+ for(i in 1:length(ns)){
+     birthday[i]<-1-exp(sum(log(1-1:(ns[i]-1)/q)))
+ }
+
+ plot(log10(ns),birthday,type="l",xlab="log10(人口)",ylab="すくなくとも同じ型のペアが1組ある確率")
+ })
   ユーザ   システム       経過  
    22.124     10.256     32.140

i = 1~1e8 まで,「1刻み」で各項に付加される確率 (q-i)/q, i=0, 1, 2, ... を計算しておき,cumprod して余事象の確率を求め,1 から引く。
さらに,この場合は対数で計算する必要はさらさらない(大きな数は出てこない)。
以下のようにすれば,計算時間は 1/3 になる

> system.time({
+ # 4.7兆分の1の確率、と考えるための数字
+ q <- 4.7e12
+
+ # 調べる人数(世界人口65億人だったり、日本人口1.2億人だったり、母集団100万人、10万人、1万人…など)
+ ns <- round(10^(seq(from=5,to=8,by=0.1)))
+ n <- 1e8
+
+ # 対数で計算する必要はない!!
+ # 同じ型のペアが少なくとも1組は存在する確率
+ birthday2 <- (1-(cumprod((q-0:(n-1))/q)))[ns]
+ plot(log10(ns), birthday2, type="l", xlab="log10(人口)", ylab="すくなくとも同じ型のペアが1組ある確率")
+ })
   ユーザ   システム       経過  
     6.135      3.136      9.195

> all.equal(birthday, birthday2)
[1] TRUE

> cbind(birthday, birthday2, birthday-birthday2)
         birthday   birthday2              
 [1,] 0.001063254 0.001063254  0.000000e+00
 [2,] 0.001684635 0.001684635  0.000000e+00
 [3,] 0.002668625 0.002668625  0.000000e+00
 [4,] 0.004226196 0.004226196  0.000000e+00
 [5,] 0.006689827 0.006689827  0.000000e+00
 [6,] 0.010581894 0.010581894  0.000000e+00
 [7,] 0.016719167 0.016719167  0.000000e+00
 [8,] 0.026368242 0.026368242 -1.110223e-16
 [9,] 0.041467409 0.041467409  0.000000e+00
[10,] 0.064919822 0.064919822  1.110223e-16
[11,] 0.100919658 0.100919658  0.000000e+00
[12,] 0.155157813 0.155157813  0.000000e+00
[13,] 0.234496703 0.234496703  0.000000e+00
[14,] 0.345260597 0.345260597  0.000000e+00
[15,] 0.488920866 0.488920866  0.000000e+00
[16,] 0.654868549 0.654868549  0.000000e+00
[17,] 0.814751459 0.814751459  0.000000e+00
[18,] 0.930901321 0.930901321  0.000000e+00
[19,] 0.985522844 0.985522844  0.000000e+00
[20,] 0.998784153 0.998784153  0.000000e+00
[21,] 0.999976020 0.999976020  0.000000e+00
[22,] 0.999999952 0.999999952  0.000000e+00
[23,] 1.000000000 1.000000000  0.000000e+00
[24,] 1.000000000 1.000000000  0.000000e+00
[25,] 1.000000000 1.000000000  0.000000e+00
[26,] 1.000000000 1.000000000  0.000000e+00
[27,] 1.000000000 1.000000000  0.000000e+00
[28,] 1.000000000 1.000000000  0.000000e+00
[29,] 1.000000000 1.000000000  0.000000e+00
[30,] 1.000000000 1.000000000  0.000000e+00
[31,] 1.000000000 1.000000000  0.000000e+00

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:もろもろ

2012年06月28日 | ブログラミング

簡易ROC」について

もろもろ。プログラムの書き方自体も相変わらずだけど(ひょっとして,Hatena では,ちゃんと書いたプログラムでも,空白がはがされてこんなになるのだろうかと??)。my.roc2 のように書く方がよいだろうということで。

my.roc<-function(x,y,t=NULL){
    Nx <- length(x)
    Ny <- length(y)
    if(is.null(t))t<-seq(from=min(x,y),to=max(x,y),length=100)
    xa<-outer(x,t,FUN="<")
    ya<-outer(y,t,FUN="<")
    xb<-apply(sign(xa),2,sum)
    yb<-apply(sign(ya),2,sum)

    xy<-data.frame(X=Nx-xb,Y=Ny-yb)
    return(xy)
}


my.roc2 <- function(x, y, t = NULL) {
    Nx <- length(x)
    Ny <- length(y)
    if (is.null(t))
        t <- seq(from = min(x, y), to = max(x, y), length = 100)
    xa <- outer(x, t, FUN = ">=") # 後で Nx-xb するなら,最初からこうしておけばよい
    ya <- outer(y, t, FUN = ">=")
    xb <- colSums(xa) # apply は遅い
    yb <- colSums(ya)
    data.frame(X = xb, Y = yb)
}

set.seed(123456)
x <- rnorm(1000, mean = 0, sd = 5)
y <- rnorm(2000, mean = 3, sd = 3)
t <- seq(from = -5, to = 20, length = 1000)
roc.out <- my.roc(x, y, t)
roc.out2 <- my.roc2(x, y, t)
all.equal(roc.out, roc.out2)
library(rbenchmark)
benchmark(my.roc(x, y, t), my.roc2(x, y, t))

結果が等しいことと,roc.out2 の方が少しは速いということで。

> all.equal(roc.out, roc.out2)
[1] TRUE
> library(rbenchmark)
> benchmark(my.roc(x, y, t), my.roc2(x, y, t))
              test replications elapsed relative user.self sys.self
1  my.roc(x, y, t)          100  42.099 1.645392    38.849    3.856
2 my.roc2(x, y, t)          100  25.586 1.000000    23.570    2.141

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

フィボナッチ数の一般項

2012年06月25日 | ブログラミング

フィボナッチ数の「一般項を求める式を用いた解法」において

当然のことではあるが,フィボナッチ数は整数のはず。示されたプログラムは,numeric(double) を使っているのだから,限界がある。しかるに,限界についてなんの言及もないのは困ったものだ。

フィボナッチ数を求めるなら,gmp ライブラリの fibnum を使うのが穏当。

> options(digits=20) # 実際は,15,6桁しか意味がない

> my.fib.3(71)
[1] 308061521170129.6875

小数点以下四捨五入だと,これ以上のフィボナッチ数は正しくない。

> fibnum(71)
Big Integer ('bigz') :
[1] 308061521170129

> my.fib.3(10)
[1] 55.000000000000014211

15,6桁しか意味がない。後ろの方の数字はゴミ。

通常,options(digits=7) だから,ゴミは丸められて表示される。

> options(digits=7)
> my.fib.3(10)
[1] 55

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:不適切な公式は使わないこと!

2012年06月15日 | ブログラミング

Rの手作りプログラム」中の「標本分散,標本共分散」について
#標本分散
sample.variance<-function(a) var(a)*(length(a)-1)/length(a)
  #以下でも同じ。
  #sample2.variance<-function(x) mean(x^2)-(mean(x))^2

#標本共分散
sample.covariance<-function(x,y) mean(x*y)-mean(x)*mean(y)

「標本分散」という使い方がおかしい。
これは詰まるところ,「変動をサンプルサイズで割る」という方の分散のこと。Excel なんかは,これは「母分散である」といっているが,これも間違い。そもそも,母分散なんか計算できるものではない。

名前は,単に「分散」という。「変動をサンプルサイズ - 1 で割る」分散は,「不偏分散」と呼ぶ。

それはさておき,mean(x^2)-(mean(x))^2 は非常にマズイ。
mean((x-mean(x))^2) とすべき。ちなみに,mean(x)^2 でよいし。
理由?それはね,場合によって計算誤差が出るから

> x <- 100000000+1:100

> sample.variance(x) # 著者が定義した関数による場合 正しい
[1] 833.25

> mean(x^2)-(mean(x))^2 # 著者が別法として提示したもの 誤り
[1] 834 # 誤った結果だ!!!

> mean((x-mean(x))^2) # 正しくは,このようにする
[1] 833.25

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:テストしてから使おう!

2012年06月15日 | ブログラミング

Rの手作りプログラム」の中の「異常値の数」
「標準偏差の±k 倍を超えた標本数」とわざわざ注記しているのに,提示されたプログラムは,その要件を満たしていない。

#異常値の数(標準偏差の±k 倍を超えた標本数)
abnormal.number <- function(x, k) {
    z <- length(x[x > k * sd(x) | x < -k * sd(x)])
    z
}

z <- length(x[x>k*sd(x) | x< -k*sd(x)]) ではなくて,
z <- length(x[x>mean(x)+k*sd(x) | x< mean(x)-k*sd(x)]) でなければならない。
さらに,これは
z <- sum(x > mean(x) + k * sd(x) | x < mean(x) - k * sd(x)) のようにすればよい(x の要素を選ぶ必要はない)。
さらに,さらに,別に結果を z に付値する必要はない。
mean(x)k*sd(x) を 2 回ずつ計算するのももったいない。
結局,プログラムは,

abnormal.number2 <- function(x, k) {
    mx <- mean(x)
    sx <- k * sd(x)
    sum(x > mx + sx | x < mx - sx)
}

ということになる。
そして,元のプログラムより若干速くなる。

> x <- rnorm(10000000)
> system.time(abnormal.number(x, 2))
   ユーザ   システム       経過  
     0.542      0.011      0.548
> system.time(abnormal.number2(x, 2))
   ユーザ   システム       経過  
     0.327      0.001      0.325

ちなみに,「標本数」というのは誤った用語法である。この場合は「サンプルサイズ」とか「標本の大きさ」の方だ。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

こんなところに計算誤差があるなんて

2012年06月15日 | ブログラミング

計算誤差については人一倍気を遣っていたつもりなのだけど,別の問題をチェックしていてキモを冷やした。

以下を見て,何とも思わない人は幸せな人なんだろうか。

> c(sin=sin(pi), cos=cos(pi), tan=tan(pi))
          sin           cos           tan
 1.224647e-16 -1.000000e+00 -1.224647e-16

> c(sin=sin(pi/2), cos=cos(pi/2), tan=tan(pi/2))
         sin          cos          tan
1.000000e+00 6.123234e-17 1.633124e+16

> c(sin=sin(pi*3/2), cos=cos(pi*3/2), tan=tan(pi*3/2))
          sin           cos           tan
-1.000000e+00 -1.836970e-16  5.443746e+15

> c(sin=sin(pi*4/2), cos=cos(pi*4/2), tan=tan(pi*4/2))
          sin           cos           tan
-2.449294e-16  1.000000e+00 -2.449294e-16

> c(sin=sin(pi*(4/2)), cos=cos(pi*(4/2)), tan=tan(pi*(4/2)))
          sin           cos           tan
-2.449294e-16  1.000000e+00 -2.449294e-16

> c(sin=sin(pi*2), cos=cos(pi*2), tan=tan(pi*2))
          sin           cos           tan
-2.449294e-16  1.000000e+00 -2.449294e-16

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

添え字をベクトルで指定するだけ

2012年06月12日 | ブログラミング

2012-06-12 ベクトルの要素にベクトル」で不思議がっているが??

「R では,ベクトルの添え字はベクトルで与えることができる。また,同じ要素を何回も重複して取り出すこともできるし,重複を許すからもとの要素数より多くの要素数を取り出すこともできる。」ということだが,どこが不思議なのかわからない。

> x <- 11:15
> x[3]
[1] 13
> x[c(3, 3, 1)]
[1] 13 13 11
> i <- c(2, 3, 1, 2, 4, 5, 4, 2, 2, 1)
> x[i]
 [1] 12 13 11 12 14 15 14 12 12 11

? "[" を見れば,"When indexing arrays by [ a single argument i can be a matrix with as many columns as there are dimensions of x; the result is then a vector with elements corresponding to the sets of indices in each row of i." と書いてある。

だから,行列の場合は行・列を示す2つの数値を二列に持つ行列で,行列の要素を取り出すことができるのである。

> (x <- matrix(11:22, 3, 4))
     [,1] [,2] [,3] [,4]
[1,]   11   14   17   20
[2,]   12   15   18   21
[3,]   13   16   19   22
> (ij <- matrix(c(1,3, 2,4, 3,2, 3,4), byrow=TRUE, nc=2))
     [,1] [,2]
[1,]    1    3
[2,]    2    4
[3,]    3    2
[4,]    3    4
> x[ij]
[1] 17 21 16 22

3次元配列では,行・列・層を表す3つの数値を3列に持つ行列で,配列の要素を取り出すことができるのである。

> (x <- array(11:34, dim=c(2, 3, 4)))
, , 1

     [,1] [,2] [,3]
[1,]   11   13   15
[2,]   12   14   16

, , 2

     [,1] [,2] [,3]
[1,]   17   19   21
[2,]   18   20   22

, , 3

     [,1] [,2] [,3]
[1,]   23   25   27
[2,]   24   26   28

, , 4

     [,1] [,2] [,3]
[1,]   29   31   33
[2,]   30   32   34

> (ijk <- matrix(c(1,3,2, 2,1,4, 2,2,2, 2,3,1), byrow=TRUE, ncol=3))
     [,1] [,2] [,3]
[1,]    1    3    2 # i=1, j=3, k=2 で,対応する要素x[1,3,2] の21が取り出されるということ
[2,]    2    1    4
[3,]    2    2    2
[4,]    2    3    1
> x[ijk]
[1] 21 30 20 16

全ては,help に書いてある。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:繰り返しをそのままプログラムすんな! その2

2012年06月12日 | ブログラミング

2012-06-12 RのLevy過程シミュレーション関数」で

プログラムというのは,規則性のある繰り返しをうまく書くべきものなのだから...

library(adehabitat)
# help(simm.levy)
system.time({
set.seed(411)
w <- simm.levy(1:500, mu = 1.5, burst = "mu = 1.5")
u <- simm.levy(1:500, mu = 2, burst = "mu = 2")
v <- simm.levy(1:500, mu = 2.5, burst = "mu = 2.5")
x <- simm.levy(1:500, mu = 3, burst = "mu = 3")
par(mfrow=c(2,2))
lapply(list(w,u,v,x), plot, perani=FALSE)
})
は,もし,パラメータを変えてもっとたくさんのシミュレーションをするときのことを考えれば,以下のように書く方がよいのは自明。
system.time({
set.seed(411)
par(mfrow=c(2,2))
mu <- seq(1.5, 3, by=0.5)
sapply(mu, function(m) plot(simm.levy(1:500, mu=m, burst=sprintf("mu = %f", m)), perani=FALSE))
})

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:数式は整理して!

2012年06月12日 | ブログラミング

R で計量ミクロ」にあるんだけど,よくまあ,こんなに複雑な式を入力したものだと思う。よく出てくる部分式をまとめて置き換えると,実に簡単な式になることがわかった。1/2 乗しているけど,そんなもの必要ない。だって,中身が2乗なのだから。2乗したものを計算してあとで平方根をとるなんて,???どうかしてるぜ~

ヘックマンのρの導出プログラム Written by Akihiko NODA †
Stata10ではヘックマンのρを直接推定していません(Stata10 Reference A-H pp.557を参照)。よって、論文に掲載する際にはちょっとした数値計算が必要となります。以下では、デルタ法(J.Wooldridge(2002) pp.44を参照)を用いてヘックマンのρを計算します。
x<- 2.5         # insert the "coef" of "/athrho"
y<- 1.2         # insert the "sd" of "/athrho"
rho<-(exp(x*2)-1)/(1+exp(x*2))
rhosd<-((exp(2*x)*2/(1+exp(2*x))-(exp(2*x)-1)*(exp(2*x)*2)/(1+exp(2*x))^2)*y^2*(exp(2*x)*2/(1+exp(2*x))-(exp(2*x)-1)*(exp(2*x)*2)/(1+exp(2*x))^2))^(1/2)
rho
rhosd

x <- 2.5  # insert the 'coef' of '/athrho'
y <- 1.2  # insert the 'sd' of '/athrho'
a <- exp(2 * x) - 1
b <- exp(2 * x) + 1
c <- exp(2 * x) * 2
rho <- a / b
rhosd <- (1 - rho) * y * c / b
rho
rhosd

同じ答えになるのはいうまでもない。すごくないだろ~~。

> x<- 2.5         # insert the "coef" of "/athrho"
> y<- 1.2         # insert the "sd" of "/athrho"
> rho<-(exp(x*2)-1)/(1+exp(x*2))
> rhosd<-((exp(2*x)*2/(1+exp(2*x))-(exp(2*x)-1)*(exp(2*x)*2)/(1+exp(2*x))^2)*y^2*(exp(2*x)*2/(1+exp(2*x))-(exp(2*x)-1)*(exp(2*x)*2)/(1+exp(2*x))^2))^(1/2)
> rho
[1] 0.9866143
> rhosd
[1] 0.03191067

> x <- 2.5  # insert the 'coef' of '/athrho'
> y <- 1.2  # insert the 'sd' of '/athrho'
> a <- exp(2 * x) - 1
> b <- exp(2 * x) + 1
> c <- exp(2 * x) * 2
> rho <- a / b
> rhosd <- (1 - rho) * y * c / b
> rho
[1] 0.9866143
> rhosd
[1] 0.03191067


コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:ダミー変数の作成

2012年06月12日 | ブログラミング

Rで計量ミクロ

> ダミー変数の作成

> Rでは説明変数に因子が含まれている場合、因子を自動でダミー変数に変換してくれるのでパッケージを利用するだけであればそのままでもかまいません。しかし、計量をSTATAやEviewsといった裏で何を計算しているのかわからない怪しいソフトウェアでなく、わざわざRでやろうとするみなさんの中には行列演算を行う場合も少なくないと思います。そんな場合は自分でダミー変数を作成する必要があります。ここでsexを、MaleとFemaleという因子をもつベクトルとします。

> > levels(sex)
> [1] "Female" "Male"

> そしてこのsexという因子ベクトルからfemaleというダミー変数を作りたいとします。 ダミー変数を作成する際は、ifelseという関数を使います。R使いなら、間違ってもforで回してifで条件分岐をしようなんて考えてはいけません。ifelse関数は、ifelse(test, yes, no)となっており、testがTUREならyesを、FALSEならnoとなるようなベクトルなり行列を返してくれる関数です。この場合、testにはFemaleならTRUEを、それ以外ならFALSEとなる条件式を入れ、yesに1を、noに0を入れればいいことになります。

> > female <- ifelse(sex=="Female",1,0)

著者も書いているように,たいていのソフトでは(R も例外ではなく),sex はカテゴリーデータだよ」と宣言しておくだけで(ダミー変数なんか作らなくても)ちゃんと分析してくれる。
しかし,めんどくさいのは,3カテゴリー以上の場合だよ。2カテゴリーの場合は簡単至極。

"Female" を 1 としたいならば,
ifelse(sex=="Female",1,0) などとしなくても female <- (sex == "Female")+0 とすればよいし,
そもそも,sex が factor なのだから,female <- 2-as.integer(sex) でよい。

実地検証

> library(rbenchmark)
> test1 <- function(x) ifelse(x == "Female", 1, 0)
> test2 <- function(x) (x == "Female")+0
> test3 <- function(x) 2-as.integer(x)
> n <- 1000000
> sex <- sample(factor(c("Male", "Female")), n, replace=TRUE)

> benchmark(test1(sex), test2(sex), test3(sex), columns = c("test", "replications", "elapsed", "relative"), replications=100)

        test replications elapsed  relative
1 test1(sex)          100  78.399 51.680290
2 test2(sex)          100  13.794  9.092947
3 test3(sex)          100   1.517  1.000000 # 50 倍速いぜ~,わいるどだろ~

> system.time(ifelse(sex == "Female", 1, 0))

   ユーザ   システム       経過  
     0.557      0.078      0.628
> system.time((sex == "Female")+0)
   ユーザ   システム       経過  
     0.115      0.007      0.122
> system.time(2-as.integer(sex))
   ユーザ   システム       経過  
      0.01       0.00       0.01

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ちょっと前の記事の数値が妥当だというシミュレーション

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

「前の記事」で指摘したことについて,その正しさを検証しておこう。

簡単なことだ。一行でできる。

> mean(replicate(1000000, sum(runif(100) < 0.1518244) >= 10))
[1] 0.950551

解説しておこうか?

runif(100) < 0.1518244

は,100 回試行して,設定した確率 0.1518244 を満たすかどうかの論理値ベクトル

sum(runif(100) < 0.1518244)

は,TRUE になった回数を数えれば,それは,アイテムを得られる回数

sum(runif(100) < 0.1518244) >= 10

は,得られたアイテムが 10 個以上であるかどうかの論理ベクトル

replicate(1000000, sum(runif(100) < 0.1518244) >= 10)

は,それを 1000000 回繰り返す。

mean(replicate(1000000, sum(runif(100) < 0.1518244) >= 10))

は,そのうち,成功した割合(0/1ベクトルの平均は,1である割合のこと)を表す。

結局,確率が 0.1518244 のときは,10 個以上のアイテムを得られる確率が 0.950551 ということ。

この数値は,試行ごとに変わるが 0.95 前後の値になることは簡単に確かめることができる。
正規近似は,近似に過ぎないこともわかる。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ブログ記事の免責について--書くのを控えていたけど

2012年06月08日 | 雑感

ガチャコンプ問題で,SNSゲームの実体というか問題点が白日の下にさらされて,それに,しぶしぶ?対応したのかとは思うが,その背景が,こんないい加減な,間違った意思決定によるものだとしたら,腹が立つ以前に笑っちゃうね。(今回は,消費者には不利でない方向で間違えているからまだまし?だけど,そう言う問題じゃない。企業側から見れば企業利益を毀損する方向の間違いだ。)

ブログの著者の所属を明らかにしてブログを書くというのは責任の所在を明らかにすると言う点では正しい選択だろうが,その記事の内容ががいい加減なもの・間違ったものだとすれば所属企業のイメージを落とすことになるよね。

自信を持って書いた記事なのだろうけど,企業はその記事の正確性について疑義があれば,著者を糾弾すべきではないかな。

よく,「この記事は,所属団体の意見とは無関係です」なんて言い逃れのおまじないを書いておくということがあるけど,この記事はそのような免責事項についての言及もないなあ。馘にされても文句は言えないのではないかな?

まあ,あんなくだらないゲームにうつつを抜かす奴の気もしれないけど,それにつけ込む奴らは許せない。

「100匹モンスターを倒せば95%の人はアイテムを10個以上手に入れられるようにする」って書いてあるけど,裏返せば「5%の人は100匹モンスターを倒してもアイテムを10個以上手に入れられない」ということですがな。モンスターを倒すのにお金が必要かは知らないが,もし,1匹のモンスターを倒すのに100円必要なら,1万円費やしても目的は達成できない人が5%いると言うことだよ。「ほとんどの人がアイテムを10個以上手に入れられるように」なんて,きれい事言っているけど,あなたが「ほとんどの人ではなかった場合の事」を考えていないと言うことが問題ではないのかということ。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ブートストラップ

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

どこかにあるけど,そこは別に最適化をしようといっているのではないので,明示しない。

for は遅い。replicate も for と同じくらい。やはり,colMeans は速い。
という,つまらない結果だけど,掲示しておく。

> library(rbenchmark)
> n <- 100000
> x <- c(1, 2, 3, 5, 3, 4, 4, 7, 8, 10, 1)

> prog1 <- function() {
+     x.boot <- numeric(n)
+     for(i in 1:n){
+         x.boot[i] <- mean(sample(x, 10, replace=T))
+     }
+     hist(x.boot)
+ }

> prog2 <- function() {
+     x.boot <- replicate(n, mean(sample(x, 10, replace=TRUE)))
+     hist(x.boot)
+ }

> prog3 <- function() {
+     x <- matrix(sample(x, n*10, replace=TRUE), 10)
+     x.boot <- colMeans(x)
+     hist(x.boot)
+ }

> benchmark(prog1(), prog2(), prog3(), replications=10)
     test replications elapsed relative user.self sys.self
1 prog1()           10  18.274 24.07642    18.140    0.195
2 prog2()           10  17.894 23.57576    17.792    0.147
3 prog3()           10   0.759  1.00000     0.697    0.054

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村