多倍長精度計算結果を見やすく出力

http://r-statistics-fan.hatenablog.com/entry/2014/08/09/103123
「Rで多倍長精度計算～とりあえず10000の階乗を計算してみる」であるが...

小数点以下も表示されてしまうので，整数部分だけを表示するように別解を。
先頭桁から区切っていくという仕様にした。別仕様は練習用。

library(Rmpfr)

one = mpfr(1,200000)
x = 10000
x = 100
fact.x = Reduce("*", c(one, 1:x))

print.mpfr = function(mpfr, n=100, m=10, Rev=FALSE) {
  num = sub("\\..+", "", formatMpfr(mpfr))
  blk = n-nchar(num) %% n
  num = c(unlist(strsplit(num, "")), rep(" ", blk))
  num = matrix(num, m)
  num = matrix(apply(num, 2, paste, collapse=""), n %/% m)
  num = apply(num, 2, paste, collapse=" ")
  num = gsub("  +", "", num)
  invisible(sapply(num, cat, "\n"))
}

> print.mpfr(fact.x) # print(fact.x) はよいが，fact.x だけではダメ
9332621544 3944152681 6992388562 6670049071 5968264381 6214685929 6389521759 9993229915 6089414639 7615651828
6253697920 8272237582 5118521091 6864000000 0000000000 00000000

なお，factorialZ() があるので，gmp パッケージも便利

print.gmp = function(gmp, n=100, m=10, Rev=FALSE) {
  num = as.character(gmp)
  blk = n-nchar(num) %% n
  num = c(unlist(strsplit(num, "")), rep(" ", blk))
  num = matrix(num, m)
  num = matrix(apply(num, 2, paste, collapse=""), n %/% m)
  num = apply(num, 2, paste, collapse=" ")
  num = gsub("  +", "", num)
  invisible(sapply(num, cat, "\n"))
}

> library(gmp)
> a = factorialZ(x)
> print.gmp(a) # print(a) ではダメ
9332621544 3944152681 6992388562 6670049071 5968264381 6214685929 6389521759 9993229915 6089414639 7615651828
6253697920 8272237582 5118521091 6864000000 0000000000 00000000

なお，文字列に直した後，gsub によって，m 桁ごとに空白，n桁ごとには改行文字を挿入してやるというのも一興かと。

print.gmp2 = function(gmp, n=100, m=10, Rev=FALSE) {
  num = as.character(gmp)
  com = sprintf('cat(gsub("(.{%i})", "\\\\1\n", gsub("(.{%i})", "\\\\1 ",num)))', (m+1)*n%/%m, m)
  eval(parse(text=com))
}

> print.gmp2(a)
9332621544 3944152681 6992388562 6670049071 5968264381 6214685929 6389521759 9993229915 6089414639 7615651828
6253697920 8272237582 5118521091 6864000000 0000000000 00000000

2 個の時系列データの相関を考えるときは...（その2）

元データがない（？）ようなので，グラフから読み取って分析してみる。

> d = read.table("takahashi.dat", header=TRUE) # データは末尾に掲載
> plot(x ~ t, type="o", col="blue", pch=16, ylim=c(-1.5, 2), ylab="x", xlab="t")
> points(y ~ t, type="o", col="red", pch=16, yaxt="n", ylab="y")
> axis(4, at=seq(-1.5, 2, length=11), labels=seq(90, 140, by=5))
> mtext("y", 4, 1.8)

> (r_xy = cor(x, y)) # 単相関係数
[1] 0.8990529
グラフから読み取ったデータなので，ちょっと誤差がある

> (r_tx = cor(t, x)) # 時間との単相関係数
[1] 0.9478756

> (r_ty = cor(t, y)) # 時間との単相関係数
[1] 0.9603337

> (r_xy.t = (r_xy-r_tx*r_ty)/sqrt(1-r_tx^2)/sqrt(1-r_ty^2)) # 偏相関係数
[1] -0.1263188     ★1★　この結果を後で参照する
このデータでは，偏相関係数は弱いながらも負の相関ということになった

> n = length(t)
> (t.stat = abs(r_xy.t)*sqrt(n-3)/sqrt(1-r_xy.t^2)) # 母偏相関係数=0 の検定統計量
[1] 0.7745721
> (p.value = pt(t.stat, df=n-3, lower.tail=FALSE)*2) # P 値を求める（両側検定）
[1] 0.4435144
有意な偏相関係数ではないということである

なお，偏相関係数は，以下のような原理で求められる。
> resid.x = lm(x ~ t)$residuals # resid.x は，時間 t で予測できない部分（時間 t に影響されない部分）
> resid.y = lm(y ~ t)$residuals # resid.y は，時間 t で予測できない部分（時間 t に影響されない部分）
> cor(resid.x, resid.y) # resid.x と resid.y の単相関係数が，すなわち，時間 t の影響を除いた，正味の x, y の相関係数
[1] -0.1263188　　　★2★　この値が，上の ★1★ で示したのと同じになる

ちなみに，resid.x で resid.y を予測するということは，時間の影響を取り除いて x で y を予測することになる
> summary(lm(resid.y ~ resid.x))

Call:
lm(formula = resid.y ~ resid.x)

Residuals:
     Min       1Q   Median       3Q      Max
-0.29687 -0.16963 -0.03915  0.13518  0.52168

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.0000     0.0314   0.000    1.000
resid.x      -0.1010     0.1286  -0.785    0.437

Residual standard error: 0.1986 on 38 degrees of freedom
Multiple R-squared:  0.01596,    Adjusted R-squared:  -0.009939
F-statistic: 0.6162 on 1 and 38 DF,  p-value: 0.4373

resid.x と resid.y を，平均値=0，分散=1に正規化して，予測式を作ると以下のようになる
> summary(lm(scale(resid.y) ~ scale(resid.x)))

Call:
lm(formula = scale(resid.y) ~ scale(resid.x))

Residuals:
    Min      1Q  Median      3Q     Max
-1.5025 -0.8585 -0.1982  0.6842  2.6403

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)
(Intercept)     7.687e-17  1.589e-01   0.000    1.000
scale(resid.x) -1.263e-01  1.609e-01  -0.785    0.437

Residual standard error: 1.005 on 38 degrees of freedom
Multiple R-squared:  0.01596,    Adjusted R-squared:  -0.009939
F-statistic: 0.6162 on 1 and 38 DF,  p-value: 0.4373

予測式の係数 scale(resid.x) が -1.263e-01 になっている。これは，★2★ の相関係数と同じ（すなわち★1★の偏相関係数と同じ）。これは，直線回帰の回帰係数と相関係数の関係からそうなる。

以上で～～すっ。

> d
    t      x      y
1   0 -1.282 -0.828
2   1 -1.273 -1.033
3   2 -1.333 -0.894
4   3 -1.308 -0.760
5   4 -1.339 -1.237
6   5 -1.020 -1.012
7   6 -0.949 -0.952
8   7 -1.096 -0.961
9   8 -1.055 -0.837
10  9 -0.920 -0.786
11 10 -0.914 -0.630
12 11 -1.001 -0.447
13 12 -0.897 -0.309
14 13 -0.915 -0.275
15 14 -0.957 -0.277
16 15 -0.887 -0.180
17 16 -0.859 -0.056
18 17 -0.531 -0.036
19 18 -0.562 -0.244
20 19 -0.718 -0.339
21 20 -0.027 -0.257
22 21  0.678 -0.245
23 22  0.235 -0.170
24 23  0.149 -0.151
25 24  0.159  0.010
26 25  0.181  0.084
27 26  0.181  0.320
28 27  0.291  0.507
29 28  0.487  0.684
30 29  0.464  0.610
31 30  0.519  0.580
32 31  0.067  0.586
33 32  0.072  0.548
34 33  0.684  0.620
35 34  0.420  0.654
36 35  0.614  0.683
37 36  0.824  0.848
38 37  0.717  1.112
39 38  0.910  1.309
40 39  1.161  1.581

2 個の時系列データの相関を考えるときは...

http://abrahamcow.hatenablog.com/entry/2014/09/11/024924
「時系列データの相関係数はあてにならない……のか？　教えて下さい」なんだけど...

私のコメントが気に触ることが多いようなのですが（特に悪意はないつもりなんですけど，すみませんね）

私は，経済学とか時系列についてはよく知らないのですが，「これは「見せかけの相関（擬似相関；spurious correlation）」の例だ」ということならば，偏相関係数を考えればよいのではないでしょうかね？？社会学などでは当たり前のように使われていると思うのですが。

> set.seed(1)
> x = cumsum(rnorm(100))
> set.seed(2)
> y = cumsum(rnorm(100))

# x，y の単相関係数
> (r_xy = cor(x, y))
[1] -0.3813307

# 時間の変数
> t = 1:100

# t，x の単相関係数
> (r_tx = cor(t, x))
[1] 0.9237727

# t，y の単相関係数
> (r_ty = cor(t, y))
[1] -0.4157964

# t を制御（固定）したときの x，y の 偏相関係数
> (r_xy.t = (r_xy-r_tx*r_ty)/sqrt(1-r_tx^2)/sqrt(1-r_ty^2))
[1] 0.00795554

ほとんど相関がないということが分かる。
検定をする場合は，通常の単相関係数の検定の計算式で検定統計量を計算し，t 分布の自由度を「サンプルサイズ - 3」にするだけ。

行列の各行から指定した要素を取り出す

k × m 行列 x と，k 個の要素を持つベクトル y がある。

> (x = matrix(1:12, 4))
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12
> y = c(1, 2, 3, 2)

x の各行 x[i,] から，y[i] 番目の要素を取り出す，簡単な方法としてはどんなのがあるか？

> apply(cbind(y, x), 1, function(z) z[z[1]+1])
[1]  1  6 11  8

追記:

このブログは，コメントは「コメント」をクリックしないと見えないので不便なので，記事中に再録する。

追記しようとここに来たところ，x[outer(y, 1:3, "==")] というのはどうかとコメントをもらっていた。

> x[outer(y, 1:3, "==")]
[1]  1  6  8 11

しかし，コメンテーターご自身もおっしゃっているように，行順でないとちょっと，困る。

うまくない理由は，R の行列が byrow=FALSE で要素が並んでいるからなので，以下のようにすればよいということになる。

> t(t(x)[outer(1:ncol(x), y, "==")])
     [,1] [,2] [,3] [,4]
[1,]    1    6   11    8

そうそう，ここにきたのは，目覚めの直前に別解がひらめいたので，追記をしにきたのだった。

> x[cbind(1:nrow(x), y)]
[1]  1  6 11  8

が素直かな。

覆面算

                INTO
                ONTO
               CANON
              INTACT
             AMMONIA
            OMISSION
           DIACRITIC
          STATISTICS
         ASSOCIATION
        ANTIMACASSAR
       CONTORTIONIST
   NONDISCRIMINATION
 + CONTRADISTINCTION
---------------------
   MISADMINISTRATION

答は，コメントとして記載しておきます