「カウント・スリー」問題

「カウント・スリー」問題

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


自然数を 1 から順に書き並べていきます。
このとき、3 の数字が現れる回数を数えます。
 
自然数 n に対し、ちょうど n 個目の 3 の数字が現れたときに書いている数を F(n) と定義します。
 
例えば F(10)＝35 です。
下の通り、10 個目の 3 は、35 を書いているときに現れます。
 
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, …
 
同様に、F(3)＝23, F(7)＝33, F(8)＝33, F(1000)＝3081 となることが確かめられます。
（「33」には 3 が 2 回現れるため、それぞれの 3 を別々に数える点に注意してください。）
 
標準入力から、自然数 n （1 ≦ n ≦ 1012）が与えられます。
標準出力に F(n) の値を出力するプログラムを書いてください。

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

「ある数までに何個の 3 があるか」という関数 g を使って，二分法（挟み撃ち法）で解を求める。

g = function(n) {
    ans = 0
    repeat {
        keta = nchar(n) - 1
        m = n %/% (10 ^ keta)
        ans = ans + m * keta * 10 ^ (keta - 1)
        if (m == 3) {
            ans = ans + n - 3 * 10 ^ keta + 1
        } else if (m > 3) {
            ans = ans + 10 ^ keta
        }
        n = n %% (10 ^ keta)
        if (n == 0) {
            break
        }
    }
    return(ans)
}

f = function(n) {
    left = 1
    right = 1000000000000
    repeat {
        middle = (left+right)%/%2
        middle.value = g(middle)-n
        if (middle.value == 0) {
            cat(middle, "\n")
            break
        }
        left.value = g(left)-n
        if (sign(left.value * middle.value) == 1) {
            left = middle
        } else {
            right = middle
        }
    }
}

# f(scan(file("stdin", "r")))

f(5) # 31
f(61) # 300
f(231) # 636
f(45609) # 88534
f(1122334455) # 1273953249
f(991357924680) # 811439945899

パズルゲーム「2048」の組み合わせは何通り？

Julia によるプログラムを追加　2021/09/25

パズルゲーム「2048」の組み合わせは何通り？

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

設問

2048というパズルゲームがあります。(Wikipedia)

iPhoneやAndroidなどのアプリだけでなく、Web上で楽しむこともできます。

実際の動かし方は試してみるとよくわかるのですが、ここでは動かし方は問いませんので、以下の遊び方は理解しなくても問題ありません。

以下、Wikipediaによる「遊び方」からの抜粋です。
4×4のマスに数字が書かれたタイルがあり、スライドさせるとそれらはマスの端まで移動し、同時に新たなタイルが出現する。
同じ数字のタイルがぶつかると2+2=4、4+4=8というように数字が足し合わされていく。
最終的に2048のタイルができればゲームクリアだが、それ以後も続けることはできる。
また、完全にタイルが動かせない状態になるとゲームオーバーとなる。

ここでは、ゲームオーバーとなった状態、つまり隣同士（各マスの上下左右）に同じ数字が並んでおらず、すべてのマスが埋まった状態だけを考えます。
このときのスコアは、それぞれのマスを作るときに足し合わされた数の和です。

例えば、「4」のマスは「2」と「2」のマスが足しあわされたもの、「8」のマスは「4」と「4」のマスが足しあわされたものですので、「4」のマスを作ると2+2=4点、「8」のマスを作るときは「4」のマスを作り、さらに足し合わせた(2+2)+(2+2)+4+4=16点、「n」のマスを作ると n×(log2n - 1)点です。この点数をすべてのマスについて足し合わせてスコアを求めます。

※iPhoneアプリなどでは最初から「4」のマスが登場する場合があるため、冒頭のスクリーンショットでは合計が一致しませんが、今回は最初にすべて「2」のマスで始まるため、「2」のマスのみ点数に加算されないものとします。（※上記の図の場合、今回の計算方法でのスコアは52536点になります。）

標準入力からスコアが与えられたとき、ゲームオーバーとなった状態が何通りあるかを求め、標準出力に出力してください。
なお、入力されるスコアは500以下の正の整数とします。
例えば、スコアが32のとき、以下の2パターンがありますので、以下のような入出力になります。

【入出力サンプル】
標準入力
32

標準出力
2

=====

n = 400 の場合，R だと 29.445 秒，Java だと 0.814 秒，C ならば 0.878 秒ということであるが，投稿先の処理系では Java は 1 秒以上で，C だと 0.4 秒で制限時間内となった。
難易度は最上位の「★ 4 つ」になっている，プログラム的にはどうということもなく，実行速度の速い言語（処理系）に依存するという結果になった。

● R

f = function(n) {
  count = 0
  x = c(0, 4, 16, 48, 128, 320)
  for (i1 in x) {
    if (i1 > n) break
    for (i2 in x) {
      i1.2 = i1 + i2
      if (i1.2 > n) break else if (i1 == i2) next
      for (i3 in x) {
        i1.3 = i1.2 + i3
        if (i1.3 > n) break else if (i2 == i3) next
        for (i4 in x) {
          i1.4 = i1.3 + i4
          if (i1.4 > n) break else if (i3 == i4) next
          for (i5 in x) {
            i1.5 = i1.4 + i5
            if (i1.5 > n) break else if (i1 == i5) next
            for (i6 in x) {
              i1.6 = i1.5 + i6
              if (i1.6 > n) break else if (i5 == i6 || i2 == i6) next
              for (i7 in x) {
                i1.7 = i1.6 + i7
                if (i1.7 > n) break else if (i3 == i7 || i6 == i7) next
                for (i8 in x) {
                  i1.8 = i1.7 + i8
                  if (i1.8 > n) break else if (i4 == i8 || i7 == i8) next
                  for (i9 in x) {
                    i1.9 = i1.8 + i9
                    if (i1.9 > n) break else if (i5 == i9) next
                    for (i10 in x) {
                      i1.10 = i1.9 + i10
                      if (i1.10 > n) break else if (i6 == i10 || i9 == i10) next
                      for (i11 in x) {
                        i1.11 = i1.10 + i11
                        if (i1.11 > n) break else if (i7 == i11 || i10 == i11) next
                        for (i12 in x) {
                          i1.12 = i1.11 + i12
                          if (i1.12 > n) break else if (i8 == i12 || i11 == i12) next
                          for (i13 in x) {
                            i1.13 = i1.12 + i13
                            if (i1.13 > n) break else if (i9 == i13) next
                            for (i14 in x) {
                              i1.14 = i1.13 + i14
                              if (i1.14 > n) break else if (i10 == i14 || i13 == i14) next
                              for (i15 in x) {
                                i1.15 = i1.14 + i15
                                if (i1.15 > n) break else if (i11 == i15 || i14 == i15) next
                                for (i16 in x) {
                                  i1.16 = i1.15 + i16
                                  if (i1.16 > n) break else if (i12 == i16 || i15 == i16) next
                                  if (i1.16 == n) count = count + 1    
  }}}}}}}}}}}}}}}}
  cat(count)
}

system.time(f(32)) # 2
system.time(f(48)) # 16
system.time(f(56)) # 56
system.time(f(256)) # 311088, 2.318sec.
system.time(f(400)) # 3430350, 28.740sec.

========================================================================================================================

● Java

import java.util.Scanner;

public class Main {

  public static void main(String[] args) {
    int n;
    if (false) {
      Scanner cin = new Scanner(System.in);
      String line;
      line = cin.nextLine();
      n = Integer.parseInt(line);
    } else {
      n = 400;
    }
    long start = System.currentTimeMillis();
    f(n);
    long end = System.currentTimeMillis();
    System.out.println((end - start) + "ms");
  }

  static void f(int n) {
    int count = 0;
    int i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16;
    int i1to2, i1to3, i1to4, i1to5, i1to6, i1to7, i1to8, i1to9, i1to10, i1to11, i1to12, i1to13, i1to14, i1to15;
    int[] x = { 0, 4, 16, 48, 128, 320 };
    for (i1 = 0; i1 < 6; i1++) {
      if (x[i1] > n) break;
      for (i2 = 0; i2 < 6; i2++) {
        i1to2 = x[i1] + x[i2];
        if (i1to2 > n) break; else if (i1 == i2) continue;
        for (i3 = 0; i3 < 6; i3++) {
          i1to3 = i1to2 + x[i3];
          if (i1to3 > n) break; else if (i2 == i3) continue;
          for (i4 = 0; i4 < 6; i4++) {
            i1to4 = i1to3 + x[i4];
            if (i1to4 > n) break; else if (i3 == i4) continue;
            for (i5 = 0; i5 < 6; i5++) {
              i1to5 = i1to4 + x[i5];
              if (i1to5 > n) break; else if (i1 == i5) continue;
              for (i6 = 0; i6 < 6; i6++) {
                i1to6 = i1to5 + x[i6];
                if (i1to6 > n) break; else if (i5 == i6 || i2 == i6) continue;
                for (i7 = 0; i7 < 6; i7++) {
                  i1to7 = i1to6 + x[i7];
                  if (i1to7 > n) break; else if (i3 == i7 || i6 == i7) continue;
                  for (i8 = 0; i8 < 6; i8++) {
                    i1to8 = i1to7 + x[i8];
                    if (i1to8 > n) break; else if (i4 == i8 || i7 == i8) continue;
                    for (i9 = 0; i9 < 6; i9++) {
                      i1to9 = i1to8 + x[i9];
                      if (i1to9 > n) break; else if (i5 == i9) continue;
                      for (i10 = 0; i10 < 6; i10++) {
                        i1to10 = i1to9 + x[i10];
                        if (i1to10 > n) break; else if (i6 == i10 || i9 == i10) continue;
                        for (i11 = 0; i11 < 6; i11++) {
                          i1to11 = i1to10 + x[i11];
                          if (i1to11 > n) break; else if (i7 == i11 || i10 == i11) continue;
                          for (i12 = 0; i12 < 6; i12++) {
                            i1to12 = i1to11 + x[i12];
                            if (i1to12 > n) break; else if (i8 == i12 || i11 == i12) continue;
                            for (i13 = 0; i13 < 6; i13++) {
                              i1to13 = i1to12 + x[i13];
                              if (i1to13 > n) break; else if (i9 == i13) continue;
                              for (i14 = 0; i14 < 6; i14++) {
                                i1to14 = i1to13 + x[i14];
                                if (i1to14 > n) break; else if (i10 == i14 || i13 == i14) continue;
                                for (i15 = 0; i15 < 6; i15++) {
                                  i1to15 = i1to14 + x[i15];
                                  if (i1to15 > n) break; else if (i11 == i15 || i14 == i15) continue;
                                  for (i16 = 0; i16 < 6; i16++) {
                                    if (i1to15 + x[i16] > n) break; else if (i12 == i16 || i15 == i16) continue;
                                    if (i1to15 + x[i16] == n) count++;
    }}}}}}}}}}}}}}}}
    System.out.println(count);
  }
}

##############

Julia では，R に比べて 40 倍ほど速い

function f(n)
  count = 0
  x = [0, 4, 16, 48, 128, 320]
  for i1 in x
    i1 > n && break
    for i2 in x
      i1_2 = i1 + i2
      i1_2 > n && break
      i1 == i2 && continue
      for i3 in x
        i1_3 = i1_2 + i3
        i1_3 > n && break
        i2 == i3 && continue
        for i4 in x
          i1_4 = i1_3 + i4
          i1_4 > n && break
          i3 == i4 && continue
          for i5 in x
            i1_5 = i1_4 + i5
            i1_5 > n && break
            i1 == i5 && continue
            for i6 in x
              i1_6 = i1_5 + i6
              i1_6 > n && break
              (i5 == i6 || i2 == i6) && continue
              for i7 in x
                i1_7 = i1_6 + i7
                i1_7 > n && break
                (i3 == i7 || i6 == i7) && continue
                for i8 in x
                  i1_8 = i1_7 + i8
                  i1_8 > n && break
                  (i4 == i8 || i7 == i8) && continue
                  for i9 in x
                    i1_9 = i1_8 + i9
                    i1_9 > n && break
                    i5 == i9 && continue
                    for i10 in x
                      i1_10 = i1_9 + i10
                      i1_10 > n && break
                      (i6 == i10 || i9 == i10) && continue
                      for i11 in x
                        i1_11 = i1_10 + i11
                        i1_11 > n && break
                        (i7 == i11 || i10 == i11) && continue
                        for i12 in x
                          i1_12 = i1_11 + i12
                          i1_12 > n && break
                          (i8 == i12 || i11 == i12) && continue
                          for i13 in x
                            i1_13 = i1_12 + i13
                            i1_13 > n && break
                            i9 == i13 && continue
                            for i14 in x
                              i1_14 = i1_13 + i14
                              i1_14 > n && break
                              (i10 == i14 || i13 == i14) && continue
                              for i15 in x
                                i1_15 = i1_14 + i15
                                i1_15 > n && break
                                (i11 == i15 || i14 == i15) && continue
                                for i16 in x
                                  i1_16 = i1_15 + i16
                                  i1_16 > n && break
                                  (i12 == i16 || i15 == i16) && continue
                                  i1_16 == n && (count += 1)
                                end
                              end
                            end
                          end
                        end
                      end
                    end
                  end
                end
              end
            end
          end
        end
      end
    end
  end
  println(count)
end

棒の長さを最小にするモビール

棒の長さを最小にするモビール

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

設問

おしゃれなインテリアとして「モビール」があります。
Wikipedia：モビール

ここでは、同じ大きさの n 個の飾りを、糸と棒を使ってバランスを取ることを考えます。
ただし、棒の長さは整数で表され、棒の長さを整数で分割した位置にしか糸を吊るせないものとします。
一般的には一つの棒の途中など複数箇所で飾りを吊るしますが、ここでは簡単にするため「棒の両端にのみ」飾りを吊るすことにします。
また、バランスを取るときに糸や棒の重さは考えなくてよいものとします。

例えば、n = 4 のとき、以下のような吊るし方が考えられます。

このとき、左の図では棒の長さの合計が6なのに対し、右の図では9となっています。
そこで、棒の長さの合計が最小になるような吊るし方を考えます。

例えば、n = 5 のときは、以下のような吊るし方をすると棒の長さの合計が11で最小になります。

標準入力から n が与えられたとき、棒の長さの合計が最小になるような吊るし方を求め、その棒の長さの合計を標準出力に出力してください。
なお、n は最大でも500以下の正の整数とします。

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

標準出力
6

==========================================================================================

以前にあった，トーナメントの作成に関連する。

n が小さい場合の状況を考察すると，
f(n) = f(i) + f(j) + n / GCM(i, j)
　　n = i + j
　　GCM は最大公約数
で，i を1 ～ int(n/2) まで変化させて，最小値を求める。
n / GCM(i, j) は，最上位の 2 つのパーツを釣り合わせるために必要な棒の長さ。

R では，採点先のシステムが弱いようで，f(499) が制限時間（1 秒）をオーバーする。

euclid = function(m, n) {
    while ((temp = n%%m) != 0) {
        n = m
        m = temp
    }
    m
}

f = function(n) {
    F = integer(n)
    F[2] = 2
    if (n >= 3) {
        for (k in 3:n) {
           Min = Inf
           for (i in 1:(k %/% 2)) {
                Min = min(Min, F[i] + F[k-i] + k / euclid(i, k-i))
           }
           F[k] = Min
        }
    }
    cat(F[n])
}
# f(scan(file("stdin", "r")))

f(50) # 107
system.time(f(499)) # 1514, 0.151sec.

awk で書き換えれば O.K. となった。

function euclid(m, n) {
      while((temp = n % m) != 0) {
          n = m
          m = temp
      }
      return m
}

function min(x, y) {
    return x < y ? x : y
}

function f(n) {
    F[1] = 0
    F[2] = 2
    if (2 >= n) {
        cat F[n]
    }
    for (k = 3; i < n; k++) {
        Min = 1e15
        for (i = 1; i
            Min = min(Min, F[i] + F[k-i] + k / euclid(i, k-i))
        }
        F[k] = Min
    }
    print F[n]
}

{
    f($1)
}

wilcox.test と wilcox_test 前者だったら残念ですね！

これも，非公開にされたのね。リンク切れになちゃった。

確かに，不都合な事実は消去したいというのはもっともだとは思うのだけど。不十分でしたと潔く振る舞うのもブシノスガタではないかな?

http://my-notes.hatenablog.com/entry/2017/09/15/173257
| My Notes
| 統計学とかR(R言語)とかPython3とかプログラミングとかの覚え書きとか走り書きとか。 座右の銘にしたい: All work and no play makes Jack a dull boy.
| 2017-09-15
| R(R言語) ノンパラメトリック検定(独立サンプルの比較、独立した2群の中心位置の比較、Mann-Whitney (マン・ウイットニー) 検定(U検定))、wilcox.test()

で，Mann-Whitney (マン・ウイットニー) 検定(U検定) -- wilcox.test において，同順位を含むデータで

> 患者_fac <- factor(c(rep("男性患者", 20), rep("女性患者", 16)))
> 評価_vec <- c(1, 4, 2, 3, 5, 3, 2, 3, 2, 5, 6, 7, 4, 2, 5, 3, 2, 2, 2, 3,
+              5, 4, 6, 4, 3, 7, 5, 6, 1, 2, 5, 6, 7, 5, 3, 5)
> df <- data.frame(患者 = 患者_fac, 評価 = 評価_vec)
> res <- wilcox.test(df$評価 ~ df$患者, data = df)
 警告メッセージ:
 wilcox.test.default(x = c(5, 4, 6, 4, 3, 7, 5, 6, 1, 2, 5, 6,  で:
   タイがあるため、正確な p 値を計算することができません
> res

    Wilcoxon rank sum test with continuity correction

data:  df$評価 by df$患者
W = 231, p-value = 0.02253
alternative hypothesis: true location shift is not equal to 0

> res$p.value
[1] 0.0225325

「タイがあるため、正確な p 値を計算することができません」に対して，「これは無視していい」とある。

wilcox.test を紹介するページの多くで「無視してよい」とあるのが無責任。せっかくノンパラメトリック検定を採用するなら，最後までちゃんとしてほしいものだ。

正確な p 値を計算する手段はある。

ひとつは exactRankTests パッケージの wilcox.exact 関数を使うやり方であった。

> library(exactRankTests)
 Package ‘exactRankTests’ is no longer under development.
 Please consider using package ‘coin’ instead.

> wilcox.exact(評価 ~ 患者, data = df)

    Exact Wilcoxon rank sum test

data:  評価 by 患者
W = 231, p-value = 0.02071
alternative hypothesis: true mu is not equal to 0

もうひとつは，library(exactRankTests) したときのメッセージ
 Package ‘exactRankTests’ is no longer under development.
 Please consider using package ‘coin’ instead.
にしたがうことである。

> library(coin)
 要求されたパッケージ survival をロード中です
> wilcox_test(評価 ~ 患者, data=df, distribution="exact")

    Exact Wilcoxon-Mann-Whitney Test

data:  評価 by 患者 (女性患者, 男性患者)
Z = 2.2974, p-value = 0.02071
alternative hypothesis: true mu is not equal to 0

wilcox.exact も wilcox_test も同じ p 値を示す。
なお，wilcox_test で表示される Z 値は p 値の計算に使われるものではないことに注意。

> pnorm(2.2974, lower.tail=FALSE)*2
[1] 0.02159596

W = 231 という統計量 W は紛らわしいが，マン・ホイットニー検定の統計量 U の大きい方である。

> n1 = 20
> n2 = 16
> (r = rank(df$評価))
 [1]  1.5 19.5  6.5 14.0 25.5 14.0  6.5 14.0  6.5 25.5 31.5 35.0 19.5  6.5 25.5 14.0  6.5  6.5
[19]  6.5 14.0 25.5 19.5 31.5 19.5 14.0 35.0 25.5 31.5  1.5  6.5 25.5 31.5 35.0 25.5 14.0 25.5
> (W1 = sum(r[1:n1])-n1*(n1+1)/2)
[1] 89
> (W2 = sum(r[n1+1:n2])-n2*(n2+1)/2)
[1] 231
> (W = max(W1, W2))
[1] 231

普通は，マン・ホイットニーの U 検定の検定統計量は W1, W2 の小さい方とされる（大きい方でもかまわないのであるが），W1+W2 = n1*n2 なので，U = n1*n2 - W で求めることができる。

> (U = n1*n2 - W)
[1] 89

示される Z を計算するためには W の期待値 E(W) と分散 V(W) が必要である。
期待値は E(U) = n1*n2/2

> (E.W = n1*n2/2) # or (W1+W2)/2
[1] 160

同順位がある場合の分散は n = n1+n2，同順位の大きさを ti として，V(W) = n1*n2/12/(n^2-n) * sum(n^3-n-Σ(ti^3-ti)

> table(r) # 同順位の個数
r
 1.5  6.5   14 19.5 25.5 31.5   35
   2    8    7    4    8    4    3
> sum(t^3-t)
[1] 1494
> n = n1+n2
> (V.W = n1*n2/12/(n^2-n)*(n^3-n-sum(t^3-t)))
[1] 955.0476

Z 値は Z= (W-E(W)) / sqrt(V(W))

> (Z = (W-E.W) / sqrt(V.W))
[1] 2.297449

しかし，wilcox_test の両側 p 値はこのZ 値から計算されているのではない。
なぜなら，

> pnorm(Z, lower.tail=FALSE)*2
[1] 0.02159318

であるから。

pnorm(Z, lower.tail=FALSE)*2 で計算される値は，連続性の補正をしないマン・ホイットニー検定の結果で示される p 値と同じものである。wilcox.test 関数を correct=FALSE で使う。

> wilcox.test(評価 ~ 患者, data = df, correct=FALSE)

    Wilcoxon rank sum test

data:  評価 by 患者
W = 231, p-value = 0.02159
alternative hypothesis: true location shift is not equal to 0

 警告メッセージ:
 wilcox.test.default(x = c(5, 4, 6, 4, 3, 7, 5, 6, 1, 2, 5, 6,  で:
   タイがあるため、正確な p 値を計算することができません

これらとは対極的に，wilcox_test では，正確な p 値が表示される。

> (a = wilcox_test(評価 ~ 患者, data=df, distribution="exact"))

    Exact Wilcoxon-Mann-Whitney Test

data:  評価 by 患者 (女性患者, 男性患者)
Z = 2.2974, p-value = 0.02071
alternative hypothesis: true mu is not equal to 0

p.value だけを取り出すためには，以下のようにする。

> pvalue(a)
[1] 0.02071037

その計算方法は，フィッシャーの正確検定 fisher.test と同じやり方で，途中の検定統計量をカイ二乗統計量ではなく，マン・ホイットニー統計量（U または W）を使う。
関数の実装例は，以下にもある。

マン・ホイットニーの U 検定(exact test)
http://aoki2.si.gunma-u.ac.jp/R/exact-mw.html

> exact.mw(df$評価[1:20], df$評価[21:36])
U = 89, Z = 2.29745, P 値 = 0.0215932
正確な P 値 = 0.02071036886
査察した分割表の個数は 14094 個

結論

詳細は以上のようであるが，どんな場合でも，coin ライブラリの wilcox_test を使おうということだ。
wilcox.test は使わない。

R の記事も多くなったけど，問題のあるものもある

> My Notes

> 統計学とかR(R言語)とかPython3とかプログラミングとかの覚え書きとか走り書きとか。 座右の銘にしたい: All work and no play makes Jack a dull boy.

> 2017-09-14

> R(R言語) ノンパラメトリック検定(適合度検定、カイ二乗検定、基準値との比較、一様性の検定)、chisq.test()

http://my-notes.hatenablog.com/entry/2017/09/14/001050

で，観測値が理論比に従うかの検定についての記事がある。

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

コードは冗長であるが，要するに，

あるカテゴリー変数について母比率が p = c(0.4, 0.2, 0.3, 0.1) であるとき，

標本での観察値が x = c(152, 102, 58, 188) のとき，標本比率が母比率と同じか？という検定をすること。

筆者は，

> x = c(152, 102, 58, 188)
> chisq.test(x)
    Chi-squared test for given probabilities

data:  x
X-squared = 77.728, df = 3, p-value < 2.2e-16

としてしまっている。これでは，母比率が一様である，つまり  p= c(0.25, 0.25, 0.25, 0.25) の検定をしていることになる。

> chisq.test(x, p= c(0.25, 0.25, 0.25, 0.25))

    Chi-squared test for given probabilities

data:  x
X-squared = 77.728, df = 3, p-value < 2.2e-16 # X-squared が同じ値

しかし，母比率は p = c(0.4, 0.2, 0.3, 0.1) なのだから，これは間違っている。

> chisq.test(x, p= c(0.4, 0.2, 0.3, 0.1))

    Chi-squared test for given probabilities

data:  x
X-squared = 448.87, df = 3, p-value < 2.2e-16

とすべき。

まあ，all or none の結論から言えばどちらも，帰無仮説を棄却することには成るんだが，間違いは間違い。

R を自習するということなら問題ないが，間違った自習結果を公表するのは問題がある。

些細なことではあるが，大きな問題は，「間違いを指摘しようにも，コメントを受け付ける態勢にない」ということ。

天上天下唯我独尊というのはいただけないなあ。（前には，コメントを受け付けていて，コメントしたこともあると思うのだが，それが不愉快だったのでコメント拒否にしたのかな）

円テーブルで席替え

円テーブルで席替え

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

設問

先生1人と生徒 n-1 人の、合わせて n 人が丸いテーブルに座っています。
途中で席替えをして、隣の人を変えることにしました。

例えば、n = 5 のとき、最初に左図のように座っていると、右図のように席替えをすればすべての人の隣が変わることになります。
このとき、先生は動かないものとします。

標準入力から n が与えられたとき、席替えをして両隣の人が変わるような並び順は何通りあるかを求め、その数を標準出力に出力してください。

n = 5 のときは、上記以外に右図の左右を反転したパターンがありますので、標準出力に「2」を出力します。
なお、n は n ≦ 10を満たす正の整数とします。

【入出力サンプル】
標準入力
5

標準出力
2

=======================================================

シンプルに並べ替えを行い，隣が全部違うかどうか調べ，そうならばカウントする。
なお，隣が同じかどうかというのは，隣との差を取り，その絶対値が 1 でも n-1 でもないということ。

permutations = function (n) {
    z = matrix(1)
    for (i in 2:n) {
        x = cbind(z, i)
        a = c(1:i, 1:(i - 1))
        z = matrix(0, ncol = ncol(x), nrow = i * nrow(x))
        z[1:nrow(x), ] = x
        for (j in 2:i - 1) {
            z[j * nrow(x) + 1:nrow(x), ] = x[, a[1:i + j]]
        }
    }
    z
}

f = function(n) {
    count = 0
    p = cbind(1, permutations(n-1)+1) # 1 が先生
    junk = apply(p, 1, function(x) {
        y = c(x, x[1])
        z = abs(diff(y))
        if (all(z != 1, z != n-1)) {
            count+1 ->> count
        }
    })
    cat(count)
}
# f(scan(file("stdin", "r")))
f(4) # 0
f(5) # 2
f(6) # 6
f(7) # 46
f(8) # 354
system.time(f(9)) # 3106, 0.499sec.
system.time(f(10)) # 29926, 3.955sec.

シンプルなプログラムでは，ちょっと計算時間が掛かる。
そこで，オンライン整数列大辞典を参照すると，https://oeis.org/A078603 があり，これは https://oeis.org/A002816 を 2 倍したもの。

f = function(n) {
    sum = factorial(n-1)/2 + (-1)^n
    for (k in 1:(n-1)) {
        for (j in 1:min(n-k, k)) {
            sum = sum+(-1)^k*choose(k-1, j-1)*choose(n-k, j)*n/(n-k)*factorial(n-k-1)/2*2^j
        }
    }
    cat(sum*2)
}

構造体のアライメント

構造体のアライメント

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

【概要】
構造体の中に、いくつかの要素が入っています。
最後の要素の先頭アドレスが、構造体の先頭から何バイトあとにあるのかを計算して下さい。
ただし、構造体の先頭アドレスは（幸運なことに）16の倍数になっています。

【入出力】
入力は
SDIL
のようになっています。
構造体の中にある要素を示す記号が、アドレス順に区切り文字なしで並んでいます。

各記号の意味は、下表のとおりです：
記号    サイズ(バイト数)    先頭アドレスの制限
C        1               なし
S        2               偶数アドレス
I        4               4の倍数アドレス
L        8               4の倍数アドレス
D        8               8の倍数アドレス
M        16              16の倍数アドレス

例えば、CMの様になっていた場合、先頭のアドレスにはCが示す要素が入っています。
その後には 15byte のパディング（変数としては使われない詰め物）が入っていて、Mが示す要素はその後に続くことになります。

出力は、SDIL の示す構造体の先頭アドレスから 20 バイト先に末尾の要素があるので
20
のような感じです。
構造体のサイズではなく、最後の要素の先頭アドレスであることに注意して下さい。

【例】
入力       出力    メモリイメージ(.はパディング )
SDIL      20      SS......DDDDDDDDIIIILLLLLLLL
CM        16      C...............MMMMMMMMMMMMMMMM
CSILDC    24      C.SSIIIILLLLLLLLDDDDDDDDC

【補足】
    •    不正な入力に対処する必要はありません。
    •    Lのアライメントは 4 の倍数です。8 の倍数ではありません。
    •    データのリオーダーはないものとします。
    •    入力文字列の長さは、1 文字以上 16 文字以下です。
    •    実際のアライメントについてはWikipedia の記事が参考になるかもしれません。

===============================================

f = function(s) {
    item = c("C", "S", "I", "L", "D", "M")
    size = c(1, 2, 4, 8, 8, 16)
    position = c(1, 2, 4, 4, 8, 16)
    address = 0
    for (i in 1:nchar(s)) {
        what = grep(substr(s, i, i), item)
        fill = address %% position[what]
        if (fill != 0) {
            address = address + (position[what] - fill)
        }
        address = address + (i != nchar(s))*size[what]
    }
    cat(address)    
}

# f(readLines(file("stdin", "r")))

f("SDIL") # 20
f("CM") # 16
f("CSILDC") # 24

f("M") # 0
f("C") # 0
f("IMCSCS") # 38
f("CSMSISLD") # 56
f("IMCSCS") # 38
f("CSMSISLD") # 56
f("LDCDCLLMD") # 80
f("CDCSSMMDIL") # 76
f("LCSMCLIMCDL") # 80
f("MISIMLCMDLLL") # 104
f("SCSSSLCDILLCI") # 56
f("CLIMSCMDLILISI") # 100

「ディバイド・アウト」問題

「ディバイド・アウト」問題

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

設問

自然数 n と素数 p に対し、n の階乗（n！）を p でこれ以上割り切れなくなるまで繰り返し割り、
その商をさらに p で割ったときの余りを F(n, p) と定義します。
 
例えば F(12, 5)＝4 です。
12！（＝479001600）は 5 で最大 2 回割ることができます。その商は 19160064 です。
さらにこれを 5 で割った余りは 4 です。
 
同様に、F(5, 3)＝1, F(20, 7)＝5, F(1000, 31)＝11 となることが確かめられます。
 
標準入力から、自然数 n と素数 p （1 ≦ n ≦ 10^12, 1 ≦ p ≦ 10^5）が半角空白区切りで与えられます。
標準出力に F(n, p) の値を出力するプログラムを書いてください。
 
 
================================================================================


# ウィルソンの定理より
# (n %% 2*p)!≡n! mod p

1 〜 n の整数のうち p を因数としてもたない数を掛け合わせる求める関数を g(n, p) と表記すると，
n = 916554098334, p = 99961 では，
g(916554098334, 99961) = g(916554098334 %% (2*99961), 99961) = g(93858, 99961)　…… (a)

1 〜 n の整数のうち，99961 を因数として持つ数は 99961*(1, 2, ..., 916554098334 %/% 99961) = 9961*(1, 2, ..., 9169116)

n = 9169116 として，
1 〜 n について，再度 g(9169116, 99961) = g(9169116 %% (2*99961), 99961) = g(172626, 99961)　…… (b)
1 〜 n の整数のうち，99961 を因数として持つ数は 99961*(1, 2, ..., 9169116 %/% 99961) = 9961*(1, 2, ..., 91)

n = 91 として，再度 g(91, 99961) であるが，これは，99961 の因数は含まれないので，普通の階乗　…… (c)

(a), (b), (c) を掛け合わせたものが g(916554098334, 99961)。
なお，途中の掛算は p を法として行う。

f = function(n, p) {
    ans = 1
    while (n != 0) {
        for (i in 1:(n %% (2*p))) {
            if (i %% p != 0) {
                ans = (ans*i) %% p
            }
        }
        n = n %/% p
    }
    cat(ans)
}

# arg = scan(file("stdin", "r"))
# f(arg[1], arg[2])

f(5, 3) # 1
f(30, 11) # 10
f(130, 13) # 6
f(321200, 307) # 124
f(8800, 99397) # 23239
f(98765432170, 17) # 8
f(916554098334, 99961) # 18502

学力テストの主成分分析のバイプロット

奥村先生の，学力テストの主成分分析のバイプロット

https://www.scoopnest.com/ja/user/h_okumura/903961406004805632-

に対して，

Taylor さんから，

> すごい。この図を一目見ただけで、視覚的に(←重要)とても多くのことが分かるってのはまさに主成分分析の威力か。あえて言うなら、第１軸、第２軸の寄与率が多少気になるところ。

というコメントがあってですね。

寄与率は
> d = read.csv("atest2017.csv")[,-1]
> rownames(d) = c("北海道", "青森", "岩手", "宮城", "秋田",
+ "山形", "福島", "茨城", "栃木", "群馬",
+ "埼玉", "千葉", "東京", "神奈川", "新潟",
+ "富山", "石川", "福井", "山梨", "長野",
+ "岐阜", "静岡", "愛知", "三重", "滋賀",
+ "京都", "大阪", "兵庫", "奈良", "和歌山",
+ "鳥取", "島根", "岡山", "広島", "山口",
+ "徳島", "香川", "愛媛", "高知", "福岡",
+ "佐賀", "長崎", "熊本", "大分", "宮崎",
+ "鹿児島", "沖縄")
> a = prcomp(d) # scale=FALSE
> print(a)
Standard deviations (1, .., p=8):
[1] 4.8600239 2.4245367 1.4255179 0.9540044 0.6784466 0.6582192 0.4617722 0.3711827

Rotation (n x k) = (8 x 8):
           PC1        PC2        PC3        PC4         PC5        PC6         PC7         PC8
小国A 0.2685353  0.2966870  0.3536034 -0.3922446 -0.47633410  0.5021088 -0.26976989  0.09587030
小国B 0.3559633  0.3486860  0.1517456  0.6129057 -0.36436550 -0.4430136 -0.09111218  0.12531363
小算A 0.2876291  0.5301904 -0.1242260 -0.4727719  0.47286916 -0.3866795 -0.13117396  0.08300222
小算B 0.3636727  0.3189432 -0.5000904  0.2562970  0.09810079  0.5004983  0.35105253 -0.25860909
中国A 0.3181836 -0.1678676  0.3342071 -0.2460678 -0.09212964 -0.1682090  0.80864133  0.08748155
中国B 0.3775471 -0.2166936  0.5331339  0.2060371  0.47944475  0.1294554 -0.22278065 -0.43308407
中数A 0.4343985 -0.4446556 -0.4062962 -0.2579537 -0.34817753 -0.2679320 -0.22925760 -0.36930797
中数B 0.3920535 -0.3670341 -0.1594133  0.1083912  0.21801803  0.1864553 -0.14962358  0.75480735

> x = a$sdev
> x^2/sum(x^2)
[1] 0.701197578 0.174510320 0.060326614 0.027018693 0.013664536 0.012861885 0.006330224 0.004090149
> cumsum(x^2)/sum(x^2)
[1] 0.7011976 0.8757079 0.9360345 0.9630532 0.9767177 0.9895796 0.9959099 1.0000000

ということで，第1，第2主成分の寄与率は，70.1% と 17.5% で，かなり違っているということです。

奥村先生の示したバイプロットは

> biplot(a, cex=0.6) # ,  pc.biplot=FALSE;  default

です。prcomp(d, scale=FALSE) であることに注意。また，biplot(a, cex=0.6, pc.biplot=FALSE) であることにも注意。

biplot(a2, cex=0.6,  pc.biplot=TRUE) とする場合には，目盛りが変わります。

> biplot(a2, cex=0.6,  pc.biplot=TRUE)

なぜ，目盛りが変わるかというと，

pc.biplot

If true, use what Gabriel (1971) refers to as a "principal component biplot", with lambda = 1 and observations scaled up by sqrt(n) and variables scaled down by sqrt(n). Then inner products between variables approximate covariances and distances between observations approximate Mahalanobis distance.

ということなのではありますが，この図からそんな情報を読み取るという意気込みを持つ人はいないと思われ，むしろ，単純に図を見て，「誤った情報をくみ取る」ことがあるのではないかと危惧されます。たとえば，「福井と中数Aが近い！！」とか。

さて，もう一方で，主成分分析は，変数を正規化するかどうか（分散共分散行列からスタートするか，相関係数行列からスタートするか）というのがあり，奥村先生は前者を取ったわけですが，後者を取った場合にはどうなるか，

> b = prcomp(d, scale=TRUE)
> print(b)
Standard deviations (1, .., p=8):
[1] 2.3667078 1.1682792 0.7013992 0.4755625 0.3574532 0.3142276 0.2452455 0.1704118

Rotation (n x k) = (8 x 8):
           PC1        PC2         PC3        PC4          PC5         PC6        PC7         PC8
小国A 0.3452283  0.3287580  0.45731718 -0.3634967 -0.621926560  0.06238154 -0.1916828  0.06356397
小国B 0.3643624  0.2958606  0.04867583  0.6916666  0.004563749 -0.50656879 -0.1764822  0.10562112
小算A 0.3093163  0.5172346 -0.12949306 -0.3864072  0.650100244  0.07318099 -0.1923574  0.07566833
小算B 0.3560964  0.2752739 -0.54310255  0.1375482 -0.279154784  0.33005870  0.4938781 -0.22969101
中国A 0.3772570 -0.2566438  0.34553984 -0.1992751  0.207000000 -0.33832347  0.6852228  0.07655498
中国B 0.3675927 -0.2743884  0.39488637  0.3058082  0.237866718  0.53923815 -0.1760226 -0.40223179
中数A 0.3433546 -0.4017514 -0.37533376 -0.2957626 -0.101814031 -0.39864743 -0.3525097 -0.44587388
中数B 0.3608830 -0.3981434 -0.24607800  0.0230828 -0.053227183  0.24706694 -0.1625342  0.74824164
> biplot(b, cex=0.6, pc.biplot=TRUE)

あまり変わらないと言えば変わらないが。

しかし，いずれにしても，biplot は ? biplot で見るように

There are many variations on biplots (see the references) and perhaps the most widely used one is implemented by biplot.princomp.

ということで，どのような基準（やりかた）で書かれたものかを理解していないと，誤った判断をしてしまう可能性が大である。

ということで，主成分（重み）ではなく，主成分負荷量と主成分得点を別々に描く方がよいのではないかと思う次第ではある。

http://aoki2.si.gunma-u.ac.jp/R/prcomp2.html にあるような，主成分分析の結果を表すときによく使われる（？）普通の表現で，

> prcomp2(d)
                　PC1     PC2     Contribution
小国A            　0.817   0.384     0.815
小国B            　0.862   0.346     0.863
小算A            　0.732   0.604     0.901
小算B            　0.843   0.322     0.814
中国A            　0.893  -0.300     0.887
中国B            　0.870  -0.321     0.860
中数A            　0.813  -0.469     0.881
中数B            　0.854  -0.465     0.946
Eigen.values    　5.601   1.365
Proportion      　70.016  17.061
Cumulative.prop.　70.016  87.077

> x = t(b$sdev*t(b$rotation))
> #quartz("atest2017", width=5, height=5)
> plot(x[,1], x[,2], pch=19, xlim=c(-1,1), ylim=c(-1,1), asp=1, xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
> axis(1, at=-4:4*0.25, pos=0)
> axis(2, at=-4:4*0.25, pos=0)
> text(x[,1], x[,2], colnames(d), pos=c(2,4,4,1, 4,2,2,4))
> abline(h=0, v=0)
> theta = seq(0, 2*pi, length=1000)
> x = cos(theta)
> y = sin(theta)
> lines(x, y, lty=3, col="red")

これによって描かれるのは，主成分負荷量。赤の破線は主成分負荷量の二乗（寄与率）が 1 となる境界。これに近いほど，寄与率が高いということ。

主成分得点のプロットは
> plot(b$x[,1], b$x[,2], pch=19, cex=0.6, xlab="PC1", ylab="PC2")
> text(b$x[,1], b$x[,2], rownames(d), cex=0.6, pos=4)

最後の 2 つを 1 枚のグラフにまとめたのが biplot の一例。

しかし，2 枚の図を別々に見れば，変な誤解は生まれようがない。

まとめ

・prcomp で scale = FALSE / TRUE いずれにするか
・biplot で pc.biplot=FALSE / TRUE いずれにするか
・biplot をやめるかどうか

 

移動量が最小のハンカチ落とし

移動量が最小のハンカチ落とし

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

設問

n人に加えて鬼一人がハンカチ落としをしています。
ハンカチ落としでは、鬼以外の人が円になって座ったあと、鬼が円の外を走ります。
鬼が誰かの後ろでハンカチを落とすと、落とされた人は鬼が一周してくるまでの間に気付き、鬼を追いかけます。

なお、走る速さは全員が同じため、鬼が追い付かれることはないものとします。
また、落とされた人はすぐに気付いて追いかけ、落とした人はその場所まで一周回ってきて座るものとします。

このとき、鬼になった人は他の鬼と違う量だけ移動してハンカチを落とすことにします。
これを繰り返し、すべての位置に一度ずつハンカチを落とすことを考えます。
（すべての位置にハンカチが落ちた時点で終了します。）

最初に鬼はAの位置にハンカチを落とします。
例えば、4人の場合、1人分→2人分→3人分を順に移動して落とすと、すべての位置に一度ずつハンカチが落ちます。
同様に、3人の場合は1人分→4人分を順に移動、もしくは2人分→5人分を順に移動して落とす方法などが考えられます。

このとき、鬼の移動量が最小になるような移動方法を求め、その移動量の和を求めてください。
上記の通り、3人の場合は5(1+4)、4人のときは6(1+2+3)が最小となります。
標準入力から整数 n が与えられたとき、鬼の移動量が最小となる移動方法の移動量の和を標準出力に出力してください。
なお、n は9以下の正の整数が与えられるものとします。

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

標準出力
6

================================================================================

f = function(n) {
    library(e1071)
    if (n == 3) {
        p = permutations(4)
        p = unique(p[p[, 1] == 1, 1:2])
    } else {
        p = permutations(n)
        p = cbind(1, p[, 1:(n - 2)] + 1)
    }
    p = p[order(rowSums(p)), ] # 移動量の合計が小さい順に試す
    for (i in 1:nrow(p)) {
        index = 0 # ハンカチが落とされる場所（n を法として）
        position = logical(n) # ハンカチが落とされた場所の記録
        position[index + 1] = TRUE
        ok = TRUE
        for (j in p[i, ]) {
            index = (index + j)%%n
            if (position[index + 1]) {
                ok = FALSE
                break
            } else {
                position[index + 1] = TRUE
            }
        }
        if (ok) {
            cat(sum(p[i,])) # 小さい順に探しているので，見つかったら探索終了
            break
        }
    }
}

# f(scan(file("stdin", "r")))

f(3) # 5
f(4) # 6
f(5) # 12
f(6) # 15
f(7) # 23
f(8) # 28
f(9) # 38

n が偶数のときは sum(1:(n-1)) である。
n が奇数のときは (4-n+n^2) / 2 であろう。