母相関係数＝０の検定をブートストラップで

「ぶーつすとらっぷ」じゃあ，ないよ。

どこぞやらで，だれかが，R を使ってブートストラップによる母相関係数の分布を推定するというブログを書いておったが，ちっとも R らしくなかったので，書いてみた。

具体例のデータは，サンプルサイズ 20 の以下のようなもの。
x = c(7.8, 25.1, 36.1, 27.1, 24.5, 11.4, 27.6, 12.3, 5.1, 4, 26.3, 21.8, 35.4, 18.5, 31.6, 7, 10.8, 24.2, 25.4, 17.9)
y = c(29.3, 13.5, 28.9, 38.8, 30.4, 8.1, 19.7, 25.9, 2.5, 16.2, 10.3, 21.2, 27.6, 12.5, 29.1, 8, 8.6, 29.7, 25.9, 13.8)

このデータから，x, y それぞれで，重複を許して 20 個のデータを抽出して相関係数を計算するというのを多数回繰り返す。

x から重複を許して 20 個のデータをリサンプリングするのは
sample(x, replace=TRUE)

y から重複を許して 20 個のデータをリサンプリングするのは
sample(x, replace=TRUE)

その相関係数を計算するのは
cor(sample(x, replace=TRUE), sample(y, replace=TRUE))

それを，たとえば 50000 回繰り返すのは
sapply(1:50000, function(i) cor(sample(x, replace=TRUE), sample(y, replace=TRUE)))

その結果を要素数 50000 のベクトル result に格納するのは
result = sapply(1:50000, function(i) cor(sample(x, replace=TRUE),
                                         sample(y, replace=TRUE)))

そのヒストグラムは
hist(result)

リサンプリングデータにおける相関係数の絶対値が，観察されたデータでの相関係数 cor(x, y) より大きい割合は，
mean(abs(result) > cor(x, y))

sum(abs(result) > cor(x, y)) / 50000 としたいかもしれないが，それは mean だ！

リサンプリングによらず，観察されたデータでの母相関係数＝0の検定（両側）の p 値は
cor.test(x, y)$p.value

sample 関数は乱数発生を伴うので，以下のプログラムでの結果は，毎回異なるが，mean(abs(result) > cor(x, y)) と cor.test(x, y)$p.value の結果はほぼ等しくなる。

> x = c(7.8, 25.1, 36.1, 27.1, 24.5, 11.4, 27.6, 12.3, 5.1, 4, 26.3, 21.8, 35.4, 18.5, 31.6, 7, 10.8, 24.2, 25.4, 17.9)
> length(sample(x, replace=TRUE))
[1] 20
> length(sample(x, replace=TRUE))
[1] 20
> x = c(7.8, 25.1, 36.1, 27.1, 24.5, 11.4, 27.6, 12.3, 5.1, 4, 26.3, 21.8, 35.4, 18.5, 31.6, 7, 10.8, 24.2, 25.4, 17.9)
> y = c(29.3, 13.5, 28.9, 38.8, 30.4, 8.1, 19.7, 25.9, 2.5, 16.2, 10.3, 21.2, 27.6, 12.5, 29.1, 8, 8.6, 29.7, 25.9, 13.8)
> result = sapply(1:50000, function(i) cor(sample(x, replace=TRUE), sample(y, replace=TRUE)))
> hist(result)
> mean(abs(result) > cor(x, y))
[1] 0.00994
> cor.test(x, y)$p.value
[1] 0.009260503

並べ替えの繰り返し2

並べ替えの繰り返し2

締め切りが 2017/02/14 10:00 AM なので，その 1 分後に投稿されるように予約

設問

1〜m までの数字が1つずつ書かれた m 枚のカードがあります。
カードは一列に並べられており、一番左のカードに書かれている数字を見て、その数の枚数だけカードを左から取り、並びを逆順にする作業を繰り返します。
作業は左端が1になるまで続けます。

例えば、m = 4 のとき、最初の並びが「3 4 2 1」の順になっていると、以下の手順を繰り返して3回で停止します。

3 4 2 1
↓　1回目：左端が3なので、左の3枚を反転
2 4 3 1
↓　2回目：左端が2なので、左の2枚を反転
4 2 3 1
↓　3回目：左端が4なので、4枚を反転
1 3 2 4

標準入力から m, n という2つの数がスペースで区切って与えられたとき、与えられた m 枚のカードに対して上記と同様の作業を n 回行って停止するような最初の並びが何通りあるかを求め、標準出力に出力してください。
（なお、mは最大で9とします。）

例えば、 m = 4, n = 3 のとき、上記のほかに以下のパターンがあり、合わせて4通りありますので、標準出力に「4」を出力します。

4 2 1 3
↓
3 1 2 4
↓
2 1 3 4
↓
1 2 3 4

4 1 3 2
↓
2 3 1 4
↓
3 2 1 4
↓
1 2 3 4

2 3 4 1
↓
3 2 4 1
↓
4 2 3 1
↓
1 3 2 4

【入出力サンプル】
標準入力
4 3

標準出力
4

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

R ならば短く書けるけど，実行時間が 1 秒を超える（せいぜい 4，5 秒なんだけど）
出題者は，再帰プログラムを望んでいるのだろうけど

f = function(m, n) {
    g = function(x) {
        count = 0
        while (x[1] != 1) {
            i = 1:x[1]
            x = c(rev(x[i]), x[-i])
            count = count+1
            if (count > n) break
        }
        count == n
    }
    library(e1071)
    a = permutations(m)
    z = apply(a, 1, g)
    cat(sum(z))
}

system.time(f(4, 3)) # 4
system.time(f(6, 2)) # 120
system.time(f(8, 15)) # 277
system.time(f(9, 30)) # 1
system.time(f(9, 0)) # 40320

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

Java ならば瞬殺

import java.util.Scanner;

public class sortCards {
    public static void main(String[] args) {
        Scanner cin = new Scanner(System.in);
        String line;
        String [] line2 = new String[2];
        line = cin.nextLine();
        line2 = line.split(" ");
        int m = Integer.parseInt(line2[0]);
        int n = Integer.parseInt(line2[1]);
        int sizeZ = factorial(m);
//      System.out.println(m+" "+sizeZ);
        int [][] z = new int[sizeZ+1][m+1];
        int ans = 0;
        z = permutations(m);
        for (int i = 1; sizeZ >= i; i++) {
            int cnt = 0;
            while (z[i][1] != 1) {
                int k = z[i][1];
//              System.out.println(k+" "+k/2);
                for (int j = 1; k/2 >= k; j++) {
                    int temp = z[i][j];
                    z[i][j] = z[i][k-j+1];
                    z[i][k-j+1] = temp;
                }
                cnt++;
            }
//          System.out.println(cnt);
            if (cnt == n) ans++;
        }
        System.out.println(ans);
    }
    
    static int factorial(int n) {
        int i, result = 1;
        for (i = 1; n >= i; i++) {
            result *= i;
        }
        return result;
    }
    
    static int[][] permutations(int n) {
        int sizeZ = factorial(n);
        int sizeX = sizeZ / (n - 1);
        int[][] z = new int[sizeZ+1][n+1];
        int[][] x = new int[sizeX+1][n+1];
        int nrowZ, ncolZ, nrowX, ncolX;
        int i, i2, j, j2;
        z[1][1] = 1;
        nrowZ = ncolZ = 1;
        for (i = 2; n >= i; i++) {
            for (i2 = 1; nrowZ >= i2; i2++) {
                for (j2 = 1; ncolZ >= j2; j2++) {
                    x[i2][j2] = z[i2][j2];
                }
                x[i2][ncolZ + 1] = i;
            }
            nrowX = nrowZ;
            ncolX = ncolZ + 1;
            for (j = 1; i >= j; j++) {
                for (j2 = 1; nrowX >= j2; j2++) {
                    for (i2 = 1; ncolX >= i2; i2++) {
                        z[(j - 1) * nrowX + j2][i2] = x[j2][(j + i2 - 2) % i
                                + 1];
                    }
                }
            }
            nrowZ = i * nrowX;
            ncolZ = ncolX;
        }
        
        return z;
    }

}