2012-01-13

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

作者: Maria L..Rizzo,石井一夫,村田真樹
出版社/メーカー: オーム社
発売日: 2011/05/26
メディア: 単行本（ソフトカバー）
購入: 1人クリック: 36回
この商品を含むブログ (5件) を見る

こちらの続き
実数の計算機表現
- 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
- 提供する機能はWikiで(こちら)

library(gsl)
help(gsl)

一次元の求根法
- 一次導関数の算出を必要とするもの(ニュートン法、ニュートン-ラプソン法)
- 一次導関数の算出を必要としないもの(ブレント法)
- uniroot()関数はブレント法
- 関連話題polyroot()関数
数値積分 integrate()関数
最尤推定

library(stats4)
help(mle)

一次元最適化 optimize()
二次元最適化 optim()
- 二次形式はその中で特徴的な位置を占める
線形計画法:シンプレックス法(と、ここまで来て、ここにつながる)
ついでに
- RのReserved words(こちら)

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