Rによる統計解析ハンドブック 第4章演習問題解答

This post is solutions of chapter 4 of `A Handbook of Statistical Analyses Using R, Second Edition'.

このポストはRによる統計解析ハンドブック(第2版)の第4章の解答になります。

1

Test based on chi-square can't reject independency of region and site of lesion, but Fisher's exact test rejects that.

カイ2乗値に基づく検定では地域と病変部位は独立でないと言えないが、Fisherの正確検定では棄却される。

data("orallesions", package="HSAUR2")
chisq.test(orallesions) # p-value = 0.1096, accept
fisher.test(orallesions) # p-value = 0.01904, reject

2

labels in graph warp each other, so should I change row names of the dataframe to avoid wapping in advance?

グラフのラベルが重なるのですが、予めrownames関数でデータフレームの行名を変えなければならないのでしょうか。

library("vcd")
data("orallesions", package="HSAUR2")
assoc(orallesions)
mosaic(orallesions)

3

Difference of p-value is derived from the that of mean. To check this, boxplot is helpful. By executing sum(p.diff>0) in R, it turns out that elements of p.diff is always larger than 0.

違いは平均値の差の大きさに由来する。これを確認するためには箱ひげ図が便利である。Rでsum(p.diff>0)とすることで、すべて0より大きな値をとっていることが分かる。

library("coin")

n.test <- 100 # number of test
p.diff <- matrix(nrow=n.test, ncol=4) # matrix to store the difference of p-value
colnames(p.diff) <- c("diff=0.5", "diff=1", "diff=1.5", "diff=2")
sz <- 12 # 15 is upper limit on my computer. Be sure to make sz less than 15!
for( i in 1:n.test ){
  for( j in 1:4 ){
    g1 <- rnorm(sz, mean=0, sd=1)
    g2 <- rnorm(sz, mean=j/2, sd=1)
    lev <- c(rep(1, sz), rep(2, sz))
    ddd <- c(g1, g2)
    gs <- as.data.frame(cbind(ddd, lev))
    colnames(gs) <c("value", "gr")
    gs$gr <- factor(gs$gr)
    p.diff[i, j] <pvalue(wilcox_test(value ~ gr, data=gs, distribution=exact()))
               - pvalue(independence_test(value ~ gr, data=gs, distribution=exact()))
  }
}

boxplot(p.diff)

Rのpairs関数の出力をカスタマイズする

pairs関数(対散布図とも言う)はデフォルトの設定では黒い点でプロットするだけなので、グラフにより情報を持たせるため、グラフをカスタマイズする。対散布図は行方向、列方向に同じ数の図が描かれるが、i行j列の位置にある図とj行i列の位置にあるグラフは実質同じなので、工夫したグラフを試してみようということである。

ここではおなじみのirisデータを基に、右上のグラフには相関係数の絶対値、左下のグラフにはデータを色で分けたグラフと、直線性が強い組についてのみ回帰直線を描いてみることにする。

右上、左下のグラフを制御するのはpairs関数の引数のうち、panel、upper.panel、lower.panelである。デフォルトではpanelは散布図を描くようになっており、upper/lower.panelはそれを使うようになっているので、点による散布図が書かれることになる。upper/lower.panelの引数に与える関数を書き換えることで、グラフに手を加えることができる。

upper <- function(x, y, ...){ # 右上の散布図を設定
  oldpar <- par(usr = c(0, 1, 0, 1)) # 3行下でテキストの位置を(0.5, 0.5)にするための調整
  v <- abs(cor(x, y))
  sz <- 1 + v * 2 # テキストのサイズを直線性の強さで調整
  text(0.5, 0.5, sprintf("%.3f", v), cex = sz)
  par(oldpar)
}

lower <- function(x, y, ...){ # 左下の散布図を設定
  points(x, y, pch = 21, bg = c("red", "green3", "blue")[unclass(iris$Species)])
  if( abs(cor(x, y)) > 0.8 ){ # 直線性の強いデータの組だけ回帰式を作成
    lmobj <- lm(y~x)
    abline(lmobj)
  }
}

pairs(iris[,1:4], upper.panel=upper, lower.panel=lower)

気をつける必要があるのは、lower関数の1行目で、ここでplotを用いて書こうとするとエラーになる点である。低水準作図関数を使うようにすること。また、lower関数内で用いたpoints関数のbgパラメータのように、データをそのまますべて渡せば、すべての散布図で自動的にうまく描いてくれる。1つ1つの散布図を気にしてパラメータを変更する必要はない。

上のコードを実行すると、次のような画像が得られる。

TopCoder SRM455 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文についておおまかに説明する。

ペーチャは蜘蛛が好きである。彼は蜘蛛を長方形の格子の各セルに1匹づつ置いた。ペーチャは何年も蜘蛛について研究してきたので、すべての蜘蛛の動きを予想した。各秒の初めにすべての蜘蛛がセルから隣り合うセルのどれか1つに移動する(あるいは格子から離れる)。各セルにペーチャは蜘蛛の動きの方向を書いており、それはA[]という文字列型配列で与えられる。i番目の要素のj番目の文字はNSEWのいずれかであり、それぞれ北、南、東、西方向への移動を表している。もし格子の外側に移動したら、床に落ちる。1秒後に空いているセルの数を返せ。

私の解答はこちら。

public class SpidersOnTheGrid {

 public int find(String[] A) {
  int nRow = A.length;
  int nCol = A[0].length();
  boolean[][] isFilled = new boolean[nRow][nCol]; // 蜘蛛でそのセルが埋まっているかを記録
  for( int i=0 ; i<A.length ; i++ ){
   for( int j=0 ; j<A[0].length() ; j++ ){
    int xpos = i; // 蜘蛛の位置
    int ypos = j;
    if( A[i].charAt(j) == 'N' ){ // セルの内容に応じて蜘蛛を動かす
     xpos--;
    }else if( A[i].charAt(j) == 'S' ){
     xpos++;
    }else if( A[i].charAt(j) == 'E' ){
     ypos++;
    }else{ // West
     ypos--;
    }
    if( xpos < 0 || xpos >= nRow || ypos < 0 || ypos >= nCol ) continue; // 床に落ちる
    isFilled[xpos][ypos] = true; // 動いた先の座標をtrueに変える
   }
  }
  
  int nFilled = 0;
  for( int i=0 ; i<A.length ; i++ ){
   for( int j=0 ; j<A[0].length() ; j++ ){
    if( isFilled[i][j] == true ) nFilled++;
   }
  }
  return nRow*nCol-nFilled; // 空セルは全体から埋まっているセル数を引く
 }

}

得点は183.69/250、2回目のsubmitでシステムテストクリア。breakとcontinueを間違える凡ミスをした。

TopCoder SRM456 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文についておおまかに説明する。

アップルワードというのは、A、P、L、Eの文字のみからなり、ちょうどその順になっている単語をいう。Aは1つ、Pは2つ以上、LとEは1つである。大文字か小文字かは問題でない。例えば"apPpPPlE"はアップルワードであり、"ABCDE"はそうではない。wordというアップルワードに変えなければならない文字列型変数が与えられたとき、アップルワードに変換するのに最小限変換しなければならない文字数を返せ。もし変換できないのであれば、-1を返すようにせよ。

私の解答はこちら。

public class AppleWord {

 public int minRep(String word) {
  int len = word.length();
  if( len < 5 ) return -1;
  int nReplace = 0;
  if( word.charAt(0) != 'A' && word.charAt(0) != 'a' ) nReplace++;
  if( word.charAt(len-2) != 'L' && word.charAt(len-2) != 'l' ) nReplace++;
  if( word.charAt(len-1) != 'E' && word.charAt(len-1) != 'e' ) nReplace++;
  for( int i=1 ; i<len-2 ; i++ ){
   if( word.charAt(i) != 'P' && word.charAt(i) != 'p' ) nReplace++;
  }
  return nReplace;
 }

}

得点は222.52/250、1回のsubmitでシステムテストクリア。変換できない状況は、wordの長さが4以下のときとなる。

Rによる統計解析ハンドブック 第3章演習問題解答

This post is solutions of chapter 3 of `A Handbook of Statistical Analyses Using R, Second Edition'.

このポストはRによる統計解析ハンドブック(第2版)の第3章の解答になります。

1

To decide which is superior, I chose wilcox.test as relative merit is decided by the absolute difference from 43 meters. t.test doesn't tell us which index is better estimator. conf.int in wilcox.test is needed to get estimated value which indicates the difference of rank. In below code execution printed value(=-3.1999...) shows that feet is superior than meter though null hypothesis couldn't rejected as p-value is slightly greater than 0.05.

どちらが優れているかを決めるため、wilcox.testを選んだ。というのも、優劣は43からの絶対誤差で決まるからである。t.testはどちらの指標が優れた推定量であるかを示してはくれない。wilcox.testのconf.intは順位の違いを示す推定値を得るために必要である。p値が0.05よりわずかに大きいため帰無仮説(2群は順位差がない)は棄却できないものの、実行時に表示される値はfeetがmeterよりも優れていることを示している。

data("roomwidth", package="HSAUR2")
convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
wt <- wilcox.test(I(abs(width*convert-43))~unit, data=roomwidth, conf.int=TRUE)
print(wt$estimate) # =-3.200

2

data("water", package="HSAUR2")
by(water, water$location, function(x){ return(cor(x$mortality, x$hardness))})

3

There is a mistake in the textbook written in Japanese. The numerator of formula in the question should be (njk-Ejk)/sqrt(Ejk). Otherwise it can't take a nagative value.

日本語版テキストの公式は間違いで、分子の方は(njk-Ejk)/sqrt(Ejk)が正しい。そうしないと(N(0, 1)の分布なのに)負の値がとれないよね～

data("pistonrings", package="HSAUR2")

calc.res <- function(x){
  sta.res <- chisq.test(x)$residuals
  x.exp <- outer(rowSums(x), colSums(x))/sum(x) # creates Ejk
  rr <- rowSums(x)
  cc <- colSums(x)
  nn <- sum(x)
  adj.res <- matrix(nrow=nrow(x), ncol=ncol(x))
  for( j in 1:nrow(x) ){
    for( k in 1:ncol(x) ){
      adj.res[j,k] <- (x[j,k]-x.exp[j,k])/sqrt(x.exp[j,k])/sqrt((1-rr[j]/nn)*(1-cc[k]/nn))
    }
  }
  adj.res <- as.data.frame(adj.res)
  rownames(adj.res) <- rownames(x)
  colnames(adj.res) <- colnames(x)
  return(list("sta.res"=sta.res, "adj.res"=adj.res))
}
pis.res <- calc.res(pistonrings)

4

The result of prop.test indicates that null hypothesis couldn't be rejected because p-value is greater than 0.05.

prop.testの結果はp値が0.05より大きいことから、(比率が一定という)帰無仮説が棄却できないことを示している。

data("rearrests", package="HSAUR2")
prop.test(rearrests, correct=FALSE)

TopCoder SRM454 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文についておおまかに説明する。

整数の桁の和は10進数で書かれた桁の和である。例えば1234は1+2+3+4=10になる。さて、A、B、Cという3つの整数が与えられる。AとBの間にあるXという整数で、Xの桁の和とCの桁の和の差の絶対値が最小になるものを返せ。もし複数のXがあり得るなら、その中で最小のものを返せ。

私の解答はこちら。

public class MinimalDifference {

 public int findNumber(int A, int B, int C) {
  int target = digit_sum(C);
  int diff = Integer.MAX_VALUE;
  int ans = A;
  for( int i=A ; i<=B ; i++ ){
   int currentDiff = Math.abs(target - digit_sum(i)); // 差の絶対値の計算
   if( currentDiff < diff ){ // Xの更新作業
    diff = currentDiff;
    ans = i;
   }
  }
  return ans;
 }
 private int digit_sum(int num){ // 桁の和を計算する関数
  // while文で0になるまで10で割った余りを足し、10で割ることを繰り返す方法もある
  String snum = "" + num;
  int s = 0;
  for( int i=0 ; i<snum.length() ; i++ ){
   s += snum.charAt(i) - '0';
  }
  return s;
 }
}

得点は235.47/250、1回のsubmitでシステムテストクリア。Xの更新作業では、複数の候補から最小値を返すので、<=としないように注意。

TopCoder SRM453.5 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文についておおまかに説明する。

あなたは大工さんで、家具をいくつか作りたいと思っている。不幸なことに、すべての道具を無くしてしまった。need[]という、i番目の要素は、スペースで区切られた家具を作るのに必要な道具のリストが与えられる。すべての家具を作るのに買わなければならない異なる道具の数を返せ。

私の解答はこちら。

import java.util.HashSet;

public class ToolsBox {

 public int countTools(String[] need) {
  HashSet<String> hs = new HashSet<String>();
  for( int i=0 ; i<need.length ; i++ ){
   String[] tools = need[i].split(" "); // スペース区切りのリスト
   for( int j=0 ; j<tools.length ; j++ ){
    hs.add(tools[j]); // ハッシュに放り込む
   }
  }
  return hs.size(); // ハッシュのサイズ=買わなければならない道具の数
 }

}

得点は246.90/250、1回のsubmitでシステムテストクリア。