ぱらぱらめくる『Rによる計算機統計学』第１１章　Ｒの数値解析

Ｒによる計算機統計学

• こちらの続き
• 実数の計算機表現
• 10進数、２進数のこと
```library(sfsmisc)
digitsBase(0:12, 8) #-- octal representation
empty.dimnames(digitsBase(0:33, 2)) # binary
```
```> digitsBase(0:12, 8) #-- octal representation
Class 'basedInt'(base = 8) [1:13]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,]    0    0    0    0    0    0    0    0    1     1     1     1     1
[2,]    0    1    2    3    4    5    6    7    0     1     2     3     4
> empty.dimnames(digitsBase(0:33, 2)) # binary
Class 'basedInt'(base = 2) [1:34]

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0
0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0
0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
```
• Rの定数
```.Machine
```
```> .Machine
\$double.eps
[1] 2.220446e-16

\$double.neg.eps
[1] 1.110223e-16

\$double.xmin
[1] 2.225074e-308

\$double.xmax
[1] 1.797693e+308

\$double.base
[1] 2

\$double.digits
[1] 53

\$double.rounding
[1] 5

\$double.guard
[1] 0

\$double.ulp.digits
[1] -52

\$double.neg.ulp.digits
[1] -53

\$double.exponent
[1] 11

\$double.min.exp
[1] -1022

\$double.max.exp
[1] 1024

\$integer.max
[1] 2147483647

\$sizeof.long
[1] 4

\$sizeof.longlong
[1] 8

\$sizeof.longdouble
[1] 12

\$sizeof.pointer
[1] 4
```
• ちょっと違うか完全に同じか
```help(identical)
```
```> example(identical)

idntcl> identical(1, NULL) ## FALSE -- don't try this with ==
[1] FALSE

idntcl> identical(1, 1.)   ## TRUE in R (both are stored as doubles)
[1] TRUE

idntcl> identical(1, as.integer(1)) ## FALSE, stored as different types
[1] FALSE

idntcl> x <- 1.0; y <- 0.99999999999

idntcl> ## how to test for object equality allowing for numeric fuzz :
idntcl> (E <- all.equal(x,y))
[1] TRUE

idntcl> isTRUE(E) # which is simply defined to just use
[1] TRUE

idntcl> identical(TRUE, E)
[1] TRUE

idntcl> ## If all.equal thinks the objects are different, it returns a
idntcl> ## character string, and the above expression evaluates to FALSE
idntcl>
idntcl> ## even for unusual R objects :
idntcl> identical(.GlobalEnv, environment())
[1] TRUE

idntcl> ### ------- Pickyness Flags : -----------------------------
idntcl>
idntcl> ## the infamous example:
idntcl> identical(0., -0.) # TRUE, i.e. not differentiated
[1] TRUE

idntcl> identical(0., -0., num.eq = FALSE)
[1] FALSE

idntcl> ## similar:
idntcl> identical(NaN, -NaN) # TRUE
[1] TRUE

idntcl> identical(NaN, -NaN, single.NA=FALSE) # differ on bit-level
[1] FALSE

idntcl> ## for functions:
idntcl> f <- function(x) x

idntcl> f
function(x) x

idntcl> g <- compiler::cmpfun(f)

idntcl> g
function(x) x
<bytecode: 0x04135bd4>

idntcl> identical(f, g)
[1] TRUE

idntcl> identical(f, g, ignore.bytecode=FALSE)
[1] FALSE

idntcl> ## Don't show:
idntcl> m0 <- m <- structure(cbind(I=1, a=1:3), foo = "bar", class = "matrix")

idntcl> attributes(m0) <- rev(attributes(m))

idntcl> names(attributes(m0)) # 'dim' remains first, interestingly...
[1] "dim"      "class"    "foo"      "dimnames"

idntcl> stopifnot(identical(0, -0),     !identical(0, -0, num.eq=FALSE),
idntcl+           identical(NaN, -NaN), !identical(NaN, -NaN, single.NA=FALSE),
idntcl+           identical(m, m0),     !identical(m, m0, attrib.as.set=FALSE) )

idntcl> ## End Don't show
idntcl>
idntcl>
idntcl>
```
• 計算の順序のこと
```n<-4
(gamma((n-1)/2)/sqrt(pi) * gamma((n-2)/2))
exp(lgamma((n-1)/2)- lgamma((n-2)/2))/sqrt(pi)

n<-400
(gamma((n-1)/2)/sqrt(pi) * gamma((n-2)/2))
exp(lgamma((n-1)/2)- lgamma((n-2)/2))/sqrt(pi)
```
```> n<-4
> (gamma((n-1)/2)/sqrt(pi) * gamma((n-2)/2))
[1] 0.5
> exp(lgamma((n-1)/2)- lgamma((n-2)/2))/sqrt(pi)
[1] 0.5
>
> n<-400
> (gamma((n-1)/2)/sqrt(pi) * gamma((n-2)/2))
[1] Inf
警告メッセージ：
1:  'gammafn' 中の値が範囲を超えています
2:  'gammafn' 中の値が範囲を超えています
> exp(lgamma((n-1)/2)- lgamma((n-2)/2))/sqrt(pi)
[1] 7.953876
```
• べき級数展開による関数の評価
• テイラー展開
• C,Fortranでは乗算を避けて効率化される
• Rでは項数を与えてベクトル化して計算することで効率化される
• C,Fortranではループを回すのも効果的
• Rではベクトル化して計算するのが効率的
• GNU科学計算ライブラリ(こちら)のR化ライブラリはgsl
```library(gsl)
help(gsl)
```
```library(stats4)
help(mle)
```
```The reserved words in R's parser are

if else repeat while function for in next break

TRUE FALSE NULL Inf NaN NA NA_integer_ NA_real_ NA_complex_ NA_character_

... and ..1, ..2 etc, which are used to refer to arguments passed down from an enclosing function. ```
• NA,NaN,Inf
• NaNに関する国際的な(IEEE)取り決め(こちら)
• さらなるレファレンス(What Every Computer Scientist Should Know about Floating-Point Arithmetic ACM Computing Surveys, 23(1))(こちら)
• RのInf,0のルールは普通の数学の極限を原則としているらしい(こちら)
```In R, basically all mathematical functions (including basic Arithmetic), are supposed to work properly with +/- Inf and NaN as input or output.

The basic rule should be that calls and relations with Infs really are statements with a proper mathematical limit. ```
```> 1/Inf
[1] 0
> log(-10)
[1] NaN
警告メッセージ：
In log(-10) :  計算結果が NaN になりました
> log(0)
[1] -Inf
> log(0)
[1] -Inf
> log(-Inf)
[1] NaN
警告メッセージ：
In log(-Inf) :  計算結果が NaN になりました
> exp(-Inf)
[1] 0
> 1/0
[1] Inf
> (-1)/0
[1] -Inf
> 1/(-0)
[1] -Inf
> (-1)/(-0)
[1] Inf
> Inf+(-Inf)
[1] NaN
> Inf/Inf
[1] NaN
> (Inf-1)/(Inf-1)
[1] NaN
> sin(0)
[1] 0
> cos(0)
[1] 1
> tan(1)
[1] 1.557408
> tan(0)
[1] 0
> sin(pi/2)
[1] 1
> cos(pi/2)
[1] 6.123032e-17
> tan(pi/2)
[1] 1.633178e+16
> asin(0)
[1] 0
> asin(1)
[1] 1.570796
> acos(0)
[1] 1.570796
> acos(1)
[1] 0
> atan(0)
[1] 0
> atan(Inf)
[1] 1.570796
> atan(-Inf)
[1] -1.570796
> atan(-Inf)/pi
[1] -0.5
> atan(Inf)-pi/2
[1] 0
> 1/0
[1] Inf
> 1/(0^2)
[1] Inf
> 1/(0^3)
[1] Inf
> 1/(-0)
[1] -Inf
> (1/(-0))^2
[1] Inf
> (1/(-0))^3
[1] -Inf
> 1/(-0)^3
[1] Inf
> 1/0
[1] Inf
> 1/(0^2)
[1] Inf
> 1/(0^3)
[1] Inf
> 1/(-0)
[1] -Inf
> (1/(-0))^2
[1] Inf
> (1/(-0))^3
[1] -Inf
> 1/(-0)^3
[1] Inf
> -0+0
[1] 0
> -0+0+(-0)
[1] 0
> 1/(-0) + 1/0
[1] NaN
> 0*(-0)
[1] 0
> 0^3
[1] 0
> (-0)^3
[1] 0
```