「2014-08-27 0からR」に例を取れば，
http://d.hatena.ne.jp/ryamada/20140827
k 種類（1～kの数が書かれているとしよう），各 n 枚のカードから，
set = rep(seq_len(k), n)
m 枚取り出し
sample(set, m)
総和が p 以上
sum(sample(set, m)) > p
になると「勝ち」
これを n.rep 回繰り返し
replicate(n.rep, sum(sample(set, m)) > p)
「勝ち」の回数を数える（勝ち率を計算するなら，sum を mean にするだけ）
sum(replicate(n.rep, sum(sample(set, m)) > p))
関数定義も for ループも不要
http://minato.sip21c.org/standardbodyweight.R
「Rの関数定義と計算例」において，
中澤さんが，「つい Vectorize() を使ってしまうのだが，本当はもっと正統派なやり方があるんだろうなあ。」と書いている。
正統派かどうか分からないが，以下のように，mapply() を使うのも一法かと。
Vectorize は使ったことがなかったのだけど，関数定義を見たら，なあんだ，mapply() を使うためのラッパーなんだ。
なお，計算に使う係数は，行列などにしておき，性別や年齢を添え字にして引用するようにすると，プログラムが分かりやすくなる（と思う。LibreOfficeの場合も同じと思う。）。
# Calculation of standard body weights, obesity degree (%), and BMI
# based on the method of Japanese School Health Statistics
# http://www.mext.go.jp/component/b_menu/other/__icsFiles/afieldfile/2013/03/29/1331750_3.pdf
# (C) Minato Nakazawa, Ph.D. 13 August 2014
Age = c(5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17)
# Ba, Ga
a = cbind(c(0.386, 0.461, 0.513, 0.592, 0.687, 0.752, 0.782, 0.783, 0.815, 0.832, 0.766, 0.656, 0.672),
c(0.377, 0.458, 0.508, 0.561, 0.652, 0.73, 0.803, 0.796, 0.655, 0.594, 0.56, 0.578, 0.598))
# Bb, Gb
b = cbind(c(23.699, 32.382, 38.878, 48.804, 61.39, 70.461, 75.106, 75.642, 81.348, 83.695, 70.989, 51.822, 53.642),
c(22.75, 32.079, 38.367, 45.006, 56.992, 68.091, 78.846, 76.934, 54.234, 43.264, 37.002, 39.057, 42.339))
stwt = function(age, sex, height) {
age = age-4
a[age, sex]*height-b[age, sex]
}
# stwtv = Vectorize(stwt)
# Definition of a sample data
sampledata = data.frame(
Age = c(7, 11, 17, 11, 14, 17),
Sex = c(1, 1, 1, 2, 2, 2), # 1: Boy, 2: Girl
Height = c(115, 155, 185, 150, 150, 150),
Weight = c(20, 46, 71, 50, 50, 50))
# Calculation
sampledata$StandardWt = mapply(stwt, sampledata$Age, sampledata$Sex, sampledata$Height)
sampledata$Obese = (sampledata$Weight - sampledata$StandardWt)/sampledata$StandardWt*100
sampledata$BMI = sampledata$Weight/(sampledata$Height/100)^2
print(sampledata)
ごもっとも。しかし，その陳述の妥当性は「...」の内容によりますね。
p1=0.12, p2=0.08 の片側検定なら，検出力を幾つに設定するかによって，「統計学的に有意な差である」というために必要なサンプルサイズは図のようになる。
普通は検出力を 0.8 にすることが多いので，700 人ずつということになる（そもそも p1=0.12, p2=0.08 なんだから，700 人ずつになるわけはないけど）。
http://aoki2.si.gunma-u.ac.jp/R/power_prop_test2.html によれば，
> power.prop.test3(Pc=0.08, Pt=0.12, r=2/3, sig.level=0.1, power=0.8)
Nc Nt
585.1472 877.7208
となるのかな？
この結果を逆に言えば，700 人ずつ以上のサンプルを集めてしまうと「そりゃあ，サンプルサイズが大きいからでしょ」ということ。別の言い方をすると，「700 人の時点で有意かどうかをいえば十分」ということ。ビッグデータなんか要らないということ。そもそも，統計学はデータは独立でなくてはならない。同じ人が何回もカウントされるというのは想定していない。本当に厳密な答えが必要ならば，「700 人ずつの，独立 2 標本を用意し，標本比率を観察して，検定すればよい」ということに尽きる。
ビッグデータが既にある。サンプルサイズは両方合わせて数万。だったら，検定などやらなくてよい（まあ，検定に費用が掛かるわけではないので，念のために計算はしておこう）。母比率の差がどれくらいで，その差がどれくらいの利益に相当するかだけをみればよい。
選挙の場合は最後の 1 票まで，正確な得票数を出す必要があるけど，上のような場合には逐次的に検定を行って，統計学的に差があることがわかったらさっさと調査を打ち切ればよい。ビッグデータになるまで待つ必要はない（「放っておいたら集まってしまうんですよ」というならしょうがない）。
以下のグラフの上は，時系列データの2変数を表示したものである。縦軸は両者同じにしている。ずいぶん見にくいし，2変数の関係はよく分からないだろう。
真ん中のグラフは，2番目のデータを調整して表示したものである。実は，上の図と2番目の図の2番目のデータ「何々」は実際のデータを0.8倍したものである。両者がよく一致している（同じ動きをしている）ことが分かるだろう。
実は，2つのデータは2014/08/15～2014/07/18の日経平均株価の始値と終値で，「何々」というのは0.8×終値である。終値が0.8掛けされているのだから，素直に描くと上のような図になる。実体は1/0.8倍して描いた真ん中の図になるのだから，どちらが実体を表しているか分かるだろう。
なお，散布図は通常は原点を描かずに描かれる。ただし，単なる散布図だと，時間の順序が不明になるので，線で結んだり，データ識別子を付けたりして下の図のように描けば，時間変動と二変数の関係を同時に表示することもできる。
結論としては，「折れ線グラフに原点は不要」とか「二軸グラフは不適切」ということではなく，「図表示の目的に合わせて，適切な図を描く」と言うことに尽きるのではないか？？
目的によっては（場合によっては）「折れ線グラフにも原点は必要」だし「2軸の目盛りを適切に設定して図を描くことも必要」なのではないか？？
上のグラフを描いたプログラム
d = structure(list(始値 = c(15317.15, 15284.38, 15111.76, 15164.73,
15022.64, 15063.73, 15138.72, 15260, 15506.87, 15474.65, 15511.54,
15732.78, 15616.89, 15564.79, 15426.98, 15342.46, 15350.28, 15367.16,
15296.42, 15174.08), 終値 = c(15318.34, 15314.57, 15213.63,
15161.31, 15130.52, 14778.37, 15232.37, 15159.79, 15320.31, 15474.5,
15523.11, 15620.77, 15646.23, 15618.07, 15529.4, 15457.87, 15284.42,
15328.56, 15343.28, 15215.71)), .Names = c("終値", "何々"
), class = "data.frame", row.names = c("1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19", "20"))
d$何々 = 0.8*d$終値
layout(matrix(1:3, 3))
old = par(mgp=c(1.8, 0.8, 0), mar=c(3, 3, 1, 4))
plot(d$始値, ylim=c(11000, 15700), type="o", ylab="円")
points(d$何々, pch=19, type="o")
legend("right", legend=c("始値", "何々"), pch=c(1, 19))
plot(d$始値, ylim=c(14500, 15700), type="o", ylab="円")
par(new=TRUE)
plot(d$何々, pch=19, type="o", ylab="", yaxt="n")
axis(4)
legend("bottomright", legend=c("始値", "何々"), pch=c(1, 19))
plot(何々~始値, d, xlab="円1", ylab="円2")
lines(何々~始値, d)
text(d$始値, d$何々, 1:20, pos=4, xpd=TRUE)
par(old)
layout(1)
以下の図は，1983年から1989年までの日経平均株価のグラフである。
「折れ線グラフは，変化を表すものだから，縦軸の原点は不要」というのは，いささか乱暴（目的によるのだろうが）ではないか。ちなみに，日経平均株価は間隔尺度では無く比尺度（ただし，0ということはあり得るので，純粋の比尺度では無いだろうが，株価0というのはあり得ない（あったらこわい）値なので）。
縦軸の目盛りがちゃんと描かれていれば，初期から終期までどれくらいの変化があったのかは，暗算すれば分かる（まあ，4倍くらいになったと言うことだが）。
だからといって，以下に示す2つのグラフが同じといえる訳ではないのは明らか。下の図は，「暗算しなくても目分量で4倍くらいかな」とわかる。逆に言えば，上の図は，うっかり見ていると（そんな奴はいないだろ!という突っ込みはなしで）とても4倍の上昇とは見えないかも知れない（4倍以上の上昇と見せるための悪意があるだろう）ということ。
原点がもっと下にあるようなデータなら，両者の違いが甚だしくなるのは，火を見るより明らか。
原点よりはるか上で，ほんの少しの振れ幅の傾向を見たいということ自体，意味がないことのように思える（±数パーセントの変化が±数兆円の変化になるのなら別だが）。
本質は，「変化を見る」ということは，「時間との相関を見る」ということ，単純な関数関係になることは少ないだろうが，もっとも単純な場合は，「時間に比例するか」ということ。そういうことならば，どこにあろうが，相関係数は同じなので何の問題もない。
しかし，多くの場合は相関関係（相関係数）ではなく線形関係（切片と傾きを伴う予測式）ではないか？この場合には，切片をどうするかは本質的なもの。つまり，上の図のように切片を（ほぼ）0としてしまうか，下の図のように0ではないとするか。
データが直線関係にない場合には更に複雑になるだろう。
新しい記事「3 次元空間の単位ベクトルを表す乱数（その4）」を参照
n 次元空間に拡張するには，（その2）の PearsonProblem3 が都合がよい。というか，2 箇所にある 3 を引数（dimension としよう）で指定してやるだけだ。
PearsonProblem4 = function(n = 4, dimension = 5, dR = 0.1, trial = 1e+05) {
xyz = array(runif(dimension * n * trial, min = -1, max = 1), dim = c(dimension, n, trial))
xyz = apply(xyz, 2:3, function(xyz2) xyz2/sqrt(sum(xyz2^2)))
apply(xyz, 3, function(xyz2) sqrt(sum(rowSums(xyz2)^2)))
}
ans = PearsonProblem4(n = 2, dimension = 2, trial = 1e+06)
mean(ans)
hist(ans, nclass = 50)
で，上のように，2 次元で 2 歩進む場合をやってみると，原点と到達点の距離が 2 に近いというのが意外と多い。
プログラムを間違えたかと，図に描いてみた。
PearsonProblem5 = function(n = 4, dR = 0.1, trial = 1e+05) {
dimension = 2
theta = seq(0, 2 * pi, length = 500)
x = n * cos(theta)
y = n * sin(theta)
plot(x, y, type = "l", asp = 1)
xyz = array(runif(dimension * n * trial, min = -1, max = 1), dim = c(dimension, n, trial))
xyz = apply(xyz, 2:3, function(xyz2) xyz2/sqrt(sum(xyz2^2)))
apply(xyz, 3, function(xyz2) {
lines(c(0, cumsum(xyz2[1, ])), c(0, cumsum(xyz2[2, ])), col="gray80")
points(cumsum(xyz2[1, ])[n], cumsum(xyz2[2, ])[n], col = "red", pch = 16)
})
}
PearsonProblem5(n = 2, trial = 100)
確かに，半径2の円周の近くに達する場合が結構ある。
パラドックスかなあ。
新しい記事「3 次元空間の単位ベクトルを表す乱数（その4）」を参照
昨日のピアソン問題では，極座標において角度をランダムに選択すると仮定した。そのようにすると，直交座標におけるベクトルの要素値は一様分布しないことになる。
そこで，直交座標で一様分布するように x, y, z を [-1, 1] の一様乱数から取り出し，ノルムが 1 になるように正規化するようにした。
プログラムは非常に簡単になった（オリジナルの作者のプログラムは VBA で書かれていたので判読が難しい）。
また，n 歩進んだ後の X, Y, Z 座標にも偏りはない（平均値は 0 に近い）。
PearsonProblem3 = function(n = 4, dR = 0.1, trial = 1e+05) {
xyz = array(runif(3 * n * trial, min = -1, max = 1), dim = c(3, n, trial))
xyz = apply(xyz, 2:3, function(xyz2) xyz2/sqrt(sum(xyz2^2)))
apply(xyz, 3, function(xyz2) sqrt(sum(rowSums(xyz2)^2)))
}
ans = PearsonProblem3(n = 3, trial = 1e+06)
mean(ans)
hist(ans, nclass = 50)
結果は，下図のようになった。分布の形が違う（オリジナルの作者のものとも違う）。
このシミュレーションでは，平均値は 1.624884 であった。
新しい記事「3 次元空間の単位ベクトルを表す乱数（その4）」を参照
http://sssslide.com/www.slideshare.net/T2C_/kadai9-report-forpublic
ピアソン問題（課題9）であるが...
私がやってみると，結果がだいぶ違う。
どうも，作者が作ったシミュレーションデータに問題があるのではないかと思う。
データの要点は，3次元空間におけるランダムな方向を向いている単位ベクトルの x, y, z 座標を作るところにある。
作者のデータの作成法は，以下のようにまとめられる。
1. 3 個の一様乱数 [0, 1) ξ1, ξ2, ξ3 を生成する
2. u = 2*ξ1-1, v = 2*ξ2-1, w = ξ3 とする
3. d = u^2 + v^2 + w^2 とする
4. d ≦ 1 の場合に
x = 2*u*w/d
y = 2*v*w/d
z = (w^2-u^2-v^2)/d とする
となっている。実際には 4. は，ベクトルの標準化をしているのだから，「d ≦ 1 の場合に」という条件は不要
このように作られたベクトルの座標を表す x, y, z は，歩数 n ごとに決まり，それぞれの合計が最終地点の X, Y, Z 座標になる。X, Y, Z は互いに無相関である。しかし，ベクトル X, Y の平均値はほぼ 0 であるが，Z の平均値は「ほぼ -1」になってしまう。
trial = 100000
X = Y = Z = numeric(trial)
for (i in 1:trial) {
x = y = z = 0
for (j in 1:3) {
repeat {
xi1 = runif(1)
xi2 = runif(1)
xi3 = runif(1)
u = 2 * xi1 - 1
v = 2 * xi2 - 1
w = xi3
d = u ^ 2 + v ^ 2 + w ^ 2
if (d <= 1) break
}
x = x + 2 * u * w / d
y = y + 2 * v * w / d
z = z + (w ^ 2 - u ^ 2 - v ^ 2) / d
}
X[i] = x
Y[i] = y
Z[i] = z
}
> mean(X)
[1] -8.390827e-05
> mean(Y)
[1] 0.002259335
> mean(Z)
[1] -0.9993304 # **** !!!!!!!!!!!!!
実はこの X, Y, Z は 3 次元の（連続）ランダムウォークの結果なので，X, Y, Z の平均値が 0 でないと，結果が歪んでしまう。
のではないかなと思う。
作者に質問したいけど，コメントを付ける手立てがない。
私の書いたプログラムは，3 次元角座標の角度 2 つをランダムに生成し，角座標から直交座標 x, y, z に戻す方法を使っている。
PearsonProblem = function(n = 3, dR = 0.1, trial = 1e+05) {
angle = matrix(runif(2 * n * trial, min = 0, max = 2 * pi), 2 * n, trial)
ans = apply(angle, 2, function(ang) {
ph = ang[1:n]
th = ang[-(1:n)]
X = sum(cos(ph) * cos(th))
Y = sum(cos(ph) * sin(th))
Z = sum(sin(ph))
sqrt(X^2 + Y^2 + Z^2)
})
}
ans = PearsonProblem(n = 3, trial = 1e+06)
mean(ans)
hist(ans, nclass = 50)
3 歩，歩いたときの原点からの距離の分布は，下図のようになった。平均値は 1.617797
