mikeo_410
  1. ３．R
    1. １．Rの常套句
    2. ２．いろいろ（R）
    3. ３．座標の引数（ベクトル、matrix、bind、list）
    4. ４．plot
    5. ５．RGL
    6. ６．ラグランジュの未定乗数法(Rsolnp)
    7. ７．常微分方程式の汎用解法（deSolve）

いろいろ（R）

  1. invisible()

関数の戻り値はターミナルに表示され、変数に代入するように記述されていれば表示されません。通常は、これで便利ですが、図形を描画する関数は別です。

図形を描く関数は、図を描いて終わりなので、変数に代入しません。しかし、検証や、更に加工する目的で、座標値などの大量のデータを出力するように作られています。

DrawEllipse3D() と言う関数を作って、楕円を3D空間に描きました。描画のon/offは、引数で指定し、座標値を返すことにします。単に描画が目的の場合に、

        u <- DrawEllipse3D(...)

のように書くのは冗長です。

DrawEllipse3D() の最後は、list(x=x,y=y,z=z) としてありましたが、

        invisible(list(x=x,y=y,z=z))

とすれば、代入されない場合は表示をしないようにしてくれます。

  1. WebGL

RGLで描いた図は、

htmlwidgets::saveWidget(rglwidget(width=x,height=x), "x.html", selfcontained = FALSE)

のようにして x.html に変換できます。

x.html と、x_filesフォルダができます。x_files には x.html が読み込むJavaScriptが入っています。同じファイルのようですが、.jsファイルの名前はバージョン番号を含んでいて、x.html を作成した時点が異なればファイル名が変わっています。

原則は、.html に対応して、_files もあるように考えられているわけですが、本来は1つあれば良いものです。

逆に言えば、個別の図に関することは、.htmlにあると言うことです。

ユーザが作成したHTML内に表示するには、iframeを使うものと思います。

iframeを割り当てて、いっぱいにRGLのhtmlを表示したい訳ですが、方法が分かりません。

rglwidgetでwidth、heightを指定すると、図の方の大きさが設定できます。

また、保存に先立って、RGLで表示している図を動かして、ページが開いたときの図の状態にして置くことは有効でした。

  1. 等差数列

例えば円を描くのに、

t <- seq(0,2*pi,length.out=361)

x <- cos(t); y <- sin(t)

lines(x,y)

のようにしています。冗長に感じていました。

t <- (0:360)*(pi/180)

と、言う手がありました。

  1. 配列の確保

for(i in 1:10) { u[i] <- i*10; }

のように書くなら、u[]が必要です。

止む無く、u <- rep(NA,10) のようにして確保していました。

  1. > rep(NA,10)
  2.  [1] NA NA NA NA NA NA NA NA NA NA
  3. > vector(length=10)  
  4.  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  5. > numeric(length=10)  
  6.  [1] 0 0 0 0 0 0 0 0 0 0

他にも vector、numeric が使えるようです。c(1,2,3) は vector のようですが、vector で割り当てられるのは、logi 型です。しかし、数値を入れることができ、数値を入れると型がnumericと表示されるようになりました。

基本的な型の配列は確保できますが、複数の項目が組になっているデータを入れる場所を予め確保する方法がありません。唯一あるのはデータフレームのようです。

  1. 2次元の配列のデータの並び

t <- (0:360)*(pi/180)

x <- cos(t); y <- sin(t)

plot(x,y,t=”l”)

と書けば円を描けます。

  1. > m <- cbind(x,y)
  2. > str(m)
  3.  num [1:361, 1:2] 1 1 0.999 0.999 0.998 ...
  4.  - attr(*, "dimnames")=List of 2
  5.   ..$ : NULL
  6.   ..$ : chr [1:2] "x" "y"
  7. > plot(m,t="l")
  8. > n <- rbind(x,y)
  9. > str(n)
  10.  num [1:2, 1:361] 1 0 0.9998 0.0175 0.9994 ...
  11.  - attr(*, "dimnames")=List of 2
  12.   ..$ : chr [1:2] "x" "y"
  13.   ..$ : NULL
  14. > plot(n,t="l")

cbind(x,y) は、xがm[,1]です。rbindだとxはn[1,]です。mをplotすると円、nをplotすると直線が描かれます。

  1. > u <- matrix(c(x,y),ncol=2)
  2. > str(u)
  3.  num [1:361, 1:2] 1 1 0.999 0.999 0.998 ...
  4. > plot(u,t="l")
  5. > v <- matrix(c(x,y),nrow=2)
  6. > str(v,t="l")
  7.  num [1:2, 1:361] 1 1 0.999 0.999 0.998 ...
  8. > plot(v,t="l")

uをplotすると円、vをplotすると直線が描かれます。

  1. > x <- c(1,3,2,4)
  2. > y <- c("a","c","b","d")
  3. > u <- cbind(x,y)
  4. > str(u)
  5.  chr [1:4, 1:2] "1" "3" "2" "4" "a" "c" "b" "d"
  6.  - attr(*, "dimnames")=List of 2
  7.   ..$ : NULL
  8.   ..$ : chr [1:2] "x" "y"
  9. > u[order(u[,1]),]
  10.      x   y
  11. [1,] "1" "a"
  12. [2,] "2" "b"
  13. [3,] "3" "c"
  14. [4,] "4" "d"
  15. > v <- rbind(x,y)
  16. > str(v)
  17.  chr [1:2, 1:4] "1" "a" "3" "c" "2" "b" "4" "d"
  18.  - attr(*, "dimnames")=List of 2
  19.   ..$ : chr [1:2] "x" "y"
  20.   ..$ : NULL
  21. > v[,order(v[1,])]
  22.   [,1] [,2] [,3] [,4]
  23. x "1"  "2"  "3"  "4"
  24. y "a"  "b"  "c"  "d"
  1. byrow

数列 c(1,2,3, 4,5,6) があって、意味的に x <- c(1,2,3) 、y <- c(4,5,6) と分けられるとします。

  1. > matrix(c(1, 2, 3, 4, 5, 6), nrow = 2)
  2.      [,1] [,2] [,3]
  3. [1,]    1    3    5
  4. [2,]    2    4    6
  5. > matrix(c(1, 2, 3, 4, 5, 6), nrow = 3)
  6.      [,1] [,2]
  7. [1,]    1    4
  8. [2,]    2    5
  9. [3,]    3    6
  10. > matrix(c(1, 2, 3, 4, 5, 6), ncol = 2)
  11.      [,1] [,2]
  12. [1,]    1    4
  13. [2,]    2    5
  14. [3,]    3    6
  15. > matrix(c(1, 2, 3, 4, 5, 6), ncol = 3)
  16.      [,1] [,2] [,3]
  17. [1,]    1    3    5
  18. [2,]    2    4    6

matrix()の nrow、ncol では、x、y が 行になるような配列（3行2列）しか作れません。

byrow=TRUE とすると、x、y が 列になるような配列（2行3列）しか作れません。

  1. > matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, byrow=TRUE)
  2.      [,1] [,2] [,3]
  3. [1,]    1    2    3
  4. [2,]    4    5    6
  5. > matrix(c(1, 2, 3, 4, 5, 6), nrow = 3, byrow=TRUE)
  6.      [,1] [,2]
  7. [1,]    1    2
  8. [2,]    3    4
  9. [3,]    5    6
  10. > matrix(c(1, 2, 3, 4, 5, 6), ncol = 2, byrow=TRUE)
  11.      [,1] [,2]
  12. [1,]    1    2
  13. [2,]    3    4
  14. [3,]    5    6
  15. > matrix(c(1, 2, 3, 4, 5, 6), ncol = 3, byrow=TRUE)
  16.      [,1] [,2] [,3]
  17. [1,]    1    2    3
  18. [2,]    4    5    6

結局、Rでは、[,1]がx、[,2]がyのように使うのが普通のことのようです。

多くの場合、配列はランダムにアクセスするので、メモリ上の連続性は必要ないので、どちらの並びでも間違えなければ良いことなので、「普通」にして置けば良さそうです。

  1. 剰余

C#で -7%3 は、-1 となります。Math.DivRem(-7,3) は、商に-2、剰余に-1 を返します。

R の場合、商は-3、剰余が2になります。

  1. > -7%/%3; -7%%3
  2. [1] -3
  3. [1] 2

どちらでも、3\times{}-2-1=3\times{}-3+2=-7 で、問題ありません。

しかし、わたしには、除数3に-3を掛けて一旦 -9 になる点で違和感があります。わたしは、除数と商の積の絶対値は被除数の絶対値を超えないと思い込んでいるようです。

Rの剰余は角度を0°から360°未満に丸めるのに都合が良いようです。

例えば－60°は、-60 %% 360 = 300 となります。

C# 式の場合の剰余は、-60 のままです。

  1. リストの[[]]

data.frame は、行と列が揃った矩形なので、配列と同じように2次元の配列として要素にアクセスできます。

[,1] は1列目を指し、列ベクトルが得られます。[1,] は1行目を指します。これは data.frame が返されます。同じ行の要素は異なる型が許されます。data.frame は、ベクトルを列としてバインドしたものと考えられます。

list() もベクトルを列としてバインドしますが、ベクトルは独立で、長さの制約がありません。したがって、行を選択するといった操作はないようです。

通常、リストに含めたベクトルは名前で、L$x のようにしているので特に気が付きませんでしたが、L[1] と書いて列を取るとリストになっています。元のベクトルは [[1]] のように括弧を2重にするようです。

良く理解しないまま、unlist() を沢山書いていたような気がします。

  1. ゼロ判定

0.1を {2}^{-1}+{2}^{-2}+\cdot{}\cdot{}\cdot{}+{2}^{-n} と2進数の浮動小数点数で表そうとすれば無限小数となります。

  1. > (0.1*0.1)==0.01  
  2. [1] FALSE
  3. > 0.1*0.1-0.01
  4. [1] 1.734723475976807e-18

したがって、0.1 と 0.1 の積と 0.01 の比較が FALSE になったりします。

数値計算の上では、極端に小さな差はゼロと見なすことが必要です。

  1. > .Machine$double.eps
  2. [1] 2.220446049250313e-16

極端に小さな値の目安にはイプシロンがあります。

64ビットの浮動小数点数の仮数部は52ビットです。仮数部は、1.xxx の .xxx を2進数で表し、指数部と合わせて、仮数部の値に２の指数部乗を乗じて値を表現しています。

仮数部の最上位のビットは 1/2 、2 の -1乗を表します。最下位のビットは-52乗を表します。この値はイプシロンと同じです。

  1. > 2^(-52)
  2. [1] 2.220446049250313e-16

イプシロンの値は、2以上の数値と比較しようとすれば、2つの数の小数点位置を合わせるときに無くなってしまう値です。したがって、比較対象の値域に応じて倍率を掛けて使う必要があります。

地表の距離をmmまで表すとします。GRS80の赤道半径は 6,378,137 mです。赤道の一周は 40,075,017 m と計算されます。{log}_{2}\left({40075017\times{}1000}\right)\simeq{}35.222 と計算され、2進数で表すには 36 ビット必要な値です。

数値計算上の都合で、普通、赤道半径を1として計算し、最後に赤道半径を乗じます。

単純に考えると、\frac{40075017\cdot{}1000}{6378137}\simeq{}6283.185  と言った値を \frac{1}{6378137}\simeq{}1.568e-07 の精度で計算することになります。これも、

{log}_{2}\left({\frac{40075017\cdot{}1000}{6378137}}\right)+\left|{{log}_{2}\left({\frac{1}{6378137}}\right)}\right|\simeq{}35.222

と計算され、必要なのは36ビットです。

1000倍細かく計算すると考えると、{log}_{2}1000\simeq{}9.966 なので、46ビット必要で、6ビットは無視でき、a、bの値が等しいかどうかを、

  1. if(abs(b-a) <= (.Machine$double.eps*64)){
  2. }

と判断できます。

|b-a| がゼロかどうかを、点ではなく、値の範囲で判断するものですが、|b-a| が64を超えると、小数点を揃えた結果、ゼロとの比較と差が無くなってしまうはずです。

しかし、実際には、演算は80ビット以上で行われるので大抵うまくいくと思います。

  1. ラグランジュ補間、ニュートン補間
  1. library(pracma)
  2. x <- seq(-2*pi,2*pi,length.out=10) # 10点の(x,y)
  3. y <- sin(x)
  4. plot(x,y,pch=16,col="red")
  5. xs <- seq(-2*pi,2*pi,length.out=361) # 点に
  6. m <- lagrangeInterp(x, y, xs) # ラグランジュ補間
  7. lines(xs,m,col="blue")
  8. n <- newtonInterp(x, y, xs)   # ニュートン補間
  9. str(n)
  10. lines(xs,n,col="green")
  1. ラグランジュの未定乗数法

Nonlinear optimization using augmented Lagrange method.

Rsolnpパッケージのsolnp()

  1. 全オブジェクトの消去

Rの環境ではメニューにあるので覚えていなかった。

Visual Code の環境で困った。

 rm(list=ls(all=TRUE))

  1. with()

deSolveパッケージのode()のヘルプのサンプルは関数をwith()で記述している。

with() は、最初の引数に含まれる値に付いているラベルを、2つ目のステートメント群で変数名とするものらしい。

  1. > L <- list(x=1,y=2,z=3)
  2. > V <- c(a=4,b=5)
  3. > func <- function(l,v){
  4. +     with(as.list(c(l,v)), {
  5. +         z<-x+y+z+a+b; a<-9
  6. +         return(list(z,a))
  7. +     })
  8. + }
  9. > str(func(L,V))
  10. List of 2
  11.  $ : num 15
  12.  $ : num 9
  13. > str(L)
  14. List of 3
  15.  $ x: num 1
  16.  $ y: num 2
  17.  $ z: num 3
  18. > str(V)
  19.  Named num [1:2] 4 5
  20.  - attr(*, "names")= chr [1:2] "a" "b"

題目一覧へmikeo_410@hotmail.com(updated: 2023/01/21)