「Excel 方眼紙はこんな時に役に立つ《Excel 方眼紙を究める》」だけど，マジかよ...
ここにも，お笑いのネタが埋蔵されている。
> 実はExcelは優れた文書作成ツールであり
> 「Excel方眼紙」とすることで、使いやすさはさらに高まる
> 罫線やセルの結合などの機能を使い、自由に表を作成できる
弱点もあると認めているが、だから止めようなんてつゆほども思っていない。
> 複雑な構造の文書を作成したとき、ほかの人がどう再編集したらよいかが一見しただけでは分かりにくい。同僚や取引先と共有するなら複雑な構造は控えるなど、配慮が必要な場合もある
> 入力したデータの再利用がしづらくなる
「効果量でものを言え」といってもね...
効果量大，中，小が 0.8, 0.5, 0.2 ということは，図に示すように明らかなことなんだけど，
対照群（平均値=50，標準偏差=10）のとき，処置群（平均値=58, 55, 52, 標準偏差=10）なわけで，
対照群では平均値以上のものの割合は 50% であるが，処置群では対照群の平均値より大きいものの割合はそれぞれ 78.8%, 69.1%，57.9% と確かに多いわけではあるが...
それらは，対照群に比べて実質的にどれくらい優れていると評価できるのか？
効果量が 0.2 程度の制がん剤は意味があるのか？おそらく無いであろう。
しかも，社会学や心理学のように，サンプルサイズをふんだんに確保することなんかできないのだから。
そんな場合には，サンプルサイズを考慮して評価せざるを得ないだろう。つまり，究極的には検定だ。
appendix
上の図を描いたプログラム
draw = function(mean=50, sd=10, cp=50, col=1, angle=45, density=30, pos=2) {
lo = mean-3.5*sd
hi = mean+3.5*sd
x = seq(lo, hi, length=1000)
y = dnorm(x, mean, sd)
lines(c(x, hi, lo), c(y, 0, 0), col=col)
x = seq(cp, hi, length=1000)
y = dnorm(x, mean, sd)
polygon(c(x, hi, cp), c(y, 0, 0), col=col, angle=angle, density=density)
text(mean, 0.047, sprintf("mean=%i, sd=%i", mean, sd), col=col, pos=pos, xpd=TRUE, cex=0.8)
text(mean, 0.041, sprintf("hatched=%.3f", pnorm(50, mean=mean, sd=sd, lower.tail=FALSE)), xpd=TRUE, col=col, pos=pos, cex=0.8)
}
draw2 = function(mean=58) {
plot(0, 0, xlim=c(15, 100), ylim=c(0,0.04), type="n", xlab="score", ylab="", bty="n", yaxt="n")
draw()
draw(mean=mean, col=2, angle=135, pos=4)
text(80, 0.03, sprintf("effect size=%.1f", (mean-50)/10), cex=0.8)
}
layout(matrix(1:3, 3))
old=par(mgp=c(1.8, 0.8, 0), mar=c(2, 2, 2, 1))
draw2(58)
draw2(55)
draw2(52)
par(old)
layout(1)
奥村先生 「裏RjpWikiに目をつけられたかな」とは，深読みしすぎですよ
私は，「心理学系の学会誌の方針には賛成できない」という立場でものをいっているつもりです。
検定と区間推定は等価（情報量は後者が多いが）なので，BASP が両方を禁止するのは至極当たり前（にしては，BASPのものの言い方は奥歯にものの挟まった感があると思いましたが）。
検定より区間推定をというのは，マイルドな統計学者は BASP などよりもずっと前から言っていたこと。それはBASP的なものではなく，より情報量が多いからということで。
Martin J. Gardner, Douglas G. Altman
Statistics with Confidence -- Confidence intervals and statistical guidelines --
British Medical Journal,
Universities Press(Belfast) 1989
David F. Bauer (1972)
Constructing confidence sets using rank statistics
Journal of the American Statistical Association 67, 687-690.
誤差を確率的に評価することができるようになったのが近代統計学なのだから，「それ以前に戻せ（戻るぞ）」というのは，アナクロニズム。それで事が済む分野（心理学など）は，「どうぞどうぞ，お好きなように。私たちはあなたたちと立場が違うので，我が道を行きますよ」ということ。
奥村先生より
🔖 Basic and Applied Social Psychology誌はNHSTP（統計的仮説検定）を禁止する http://www.tandfonline.com/doi/full/10.100/01973533.2015.1012991…
posted at 11:06:58
心理学ってのは，理系じゃないから。
「ジェンダーの違いで，ナントカカントカに有意な差があったとかなかったとか」なんてことは，実生活に何の関係もないたわごと。そこに，科学的な確率論を持ち込まれると，困る人が多いというだけの話かな。
心理学系学会誌では，かなり前から「p 値使うな」的なことをいうところは結構あった。
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc002/027.html
1999/06/04 (金) 22:15
> アメリカ心理学会は，投票によって統計的有意度の使用を廃止することに決定しているそうです。実際，最近の論文は信頼区間推定を用いたものが増えています。
「エクセルを使いたければ，F1キーを抜け！？」なんだけど，なんだかなあ...
文章の端々に，お笑いの種がちりばめられている
> 投資銀行では、計算ミスは絶対に許されない
> エクセルを使ってしっかり計算
> 投資銀行は計算ミスを起こさないように、十分に時間をかけて計算チェック
> ミスのない計算をしなければならないところが、投資銀行の仕事がきついと言われる理由の1つ
> Ｆ2キーですべての計算の中身をチェック
> （F2キーは）数式をチェックするのにとても便利
> 1日に何百回も［Ｆ2］を押していると、間違って隣の［Ｆ1］を押してしまう
> ［Ｆ1］キーを押すと、ヘルプ画面が出る
> この画面がうっとうしいし、いちいち画面を閉じるのも面倒
> （そこで）［Ｆ1］キーをキーボードから抜いて、ミスタッチをしてもヘルプ画面が出ないようにする
> これだけ徹底するあたりが投資銀行らしいですね。特に投資銀行の若手がやっています。みなさんも試しにやってみてください。
誰がするかい！！
せめて，見た目が悪いから，スプリングだけ抜いて，キートップは戻しておくほうがよいと思うよね。
ほんと，だせーーー
> 計算を始めるとき、 最初に押すのは[+]キー
> その理由はキーボードにあります。［＝］を入れるには、［Shift］＋［＝］キーを押すので、2 手必要です。
ちなみに，Mac の US キーボードだと，= はシフトキー抜きで押せるけどね。Mac で US キーボードにしたら？？（ひらがなが刻印されてなくてすっきりしてるしね）
http://districtdatalabs.silvrback.com/intro-to-r-for-microsoft-excel-users
であるが。ggplot2 と dplyr をお勧めなのはいいとして，
diamonds <- data.frame(ggplot2::diamonds)
diamonds$size[diamonds$carat < 0.5] <- "Small"
diamonds$size[diamonds$carat >=0.5 & diamonds$carat < 1] <- "Medium"
diamonds$size[diamonds$carat >= 1] <- "Large"
というのを例示しているが，すじが悪い。
diamonds$size2 <- cut(diamonds$carat, c(0, 0.5, 1, 99999), right=FALSE, labels=c("Small", "Medium", "Large"))
のほうを勧めるほうがよいと思われる。
ifelse もお勧めではない。
Brunner-Munzel 検定 ッテ，もはや，参照する価値もない記事です。パチもん検定です。
後半に追記あり（2021/10/13)
http://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
で紹介されているが，結構時間がかかる。
> system.time(print(okumura(x, y)))
[1] 0.008037645
ユーザ システム 経過
1898.105 10.102 1923.545
brunner.munzel.test の必要な部分だけ取り出して，2倍強の速度になった。
> system.time(print(brunner.munzel.permutation.test(x, y)$p.value))
[1] 0.008037645
ユーザ システム 経過
788.666 5.037 794.799
brunner.munzel.permutation.test = function(x, y) {
BM = function(X) {
x = xandy[X]
y = xandy[-X]
r1 = rank(x)
r2 = rank(y)
r = rank(c(x, y))
m1 = mean(r[1:n1])
m2 = mean(r[n1 + 1:n2])
v1 = sum((r[1:n1] - r1 - m1 + (n1 + 1)/2)^2)/(n1 - 1)
v2 = sum((r[n1 + 1:n2] - r2 - m2 + (n2 + 1)/2)^2)/(n2 - 1)
(m2 - m1)/sqrt(n1 * v1 + n2 * v2)
}
x = na.omit(x)
y = na.omit(y)
n1 = length(x)
n2 = length(y)
xandy = c(x, y)
structure(list(method = "Brunner-Munzel Permutation Test",
p.value = mean(abs(combn(n1+n2, n1, BM)) >= abs(BM(1:n1))),
data.name = paste(deparse(substitute(x)), "and", deparse(substitute(y)))),
class = "htest")
}
okumura = function(x, y) {
bm = brunner.munzel.test(x, y)$statistic
n1 = length(x)
n2 = length(y)
N = n1 + n2
xandy = c(x, y)
foo = function(X) {
brunner.munzel.test(xandy[X], xandy[-X])$statistic
}
z = combn(1:N, n1, foo)
mean(abs(z) >= abs(bm))
}
########################################################
Julia で書いたら，マシンも違うが，約10秒（少なくとも35倍ほど速い）。
using Combinatorics # combinations
using StatsBase # tiedrank
using Statistics # mean
function brunner_munzel_permutation_test(x, y)
function BM(n1, n2, n, x, y)
r1 = tiedrank(x)
r2 = tiedrank(y)
r = tiedrank(vcat(x, y))
m1 = mean(r[1:n1])
m2 = mean(r[n1 + 1:n])
v1 = sum((r[1:n1] .- r1 .- m1 .+ (n1 + 1)/2).^2)/(n1 - 1)
v2 = sum((r[n1 + 1:n] .- r2 .- m2 .+ (n2 + 1)/2).^2)/(n2 - 1)
abs(m2 - m1)/sqrt(n1 * v1 + n2 * v2)
end
all_perm(x, n) = Iterators.product((x for i = 1:n)...)
x = [z for z in x if !ismissing(z)]; # 欠損値 missing を除くデータ
y = [z for z in y if !ismissing(z)]; # 欠損値 missing を除くデータ
n1 = length(x)
n2 = length(y)
xandy = vcat(x, y);
n = n1 + n2
s0 = BM(n1, n2, n, x, y)
results = Float64[]
for i in combinations(1:n, n1)
flag = falses(n)
flag[i] .= true
x = xandy[flag]
y = xandy[.!flag]
append!(results, BM(n1, n2, n, x, y))
end
mean(results .>= s0 - eps())
end
暗号を解けと...
「原宿のここに今度の日曜の10時だって」
差し出されたスマートフォンの画面には
69, 103, 103, 115, 39, 110, 32, 84, 104, 105, 110, 103, 115
という数字の列があった。
「ごめんね、なんかめんどくさくって。」
「麻衣の周りには不思議と暗号好きが多いねｗ」
「うんホントだね、なんでだろう。でも、麻衣はぜんぜんめんどくさくないから嫌いにならないでねっ」
「うーん、これだけじゃわからないな。ハルはなにかヒント言ってなかった？」
「ああ、ごめん。『ヒントは"MONSTER SURPRISED YOU"。 ひらがなじゃないけど日本語版だからね。』って言ってた。」
ええ？ヒントは目くらましか？？
というか，ヒントなしで 1 秒で分かった。へへへ。
intToUtf8(c(69, 103, 103, 115, 39, 110, 32, 84, 104, 105, 110, 103, 115))
"Eggs'n Things" となるがねえ http://www.eggsnthingsjapan.com/
"MONSTER SURPRISED YOU" の意味が分かんね～
レールフェンス暗号は入力文字の1文字目から1文字とばしながら上段にあたる文字を取り出して暗号文の前半を作ったあと、入力文字の2文字目から1文字とばしながら下段にあたる文字を取り出して暗号文の後半を作る暗号のこと。
以下のレールフェンス暗号を平文に戻せ
hnggdtraaeaeeu
わお。恥ずかしい。
func = function(s) {
paste(c(matrix(unlist(strsplit(s, "")), 2, byrow=TRUE)), collapse="")
}
func("hnggdtraaeaeeu")
2015/02/09 19:11 NHK の RDD による世論調査
調査対象 1496 人，回答者 978 人。つまり，回答率 0.65（低い）
安倍内閣を「支持する」と答えた人は、先月より 4 ポイント上がって 54% でした。
今月は，支持する:0.54*978 ≒ 528 人，支持しない:978-528 = 450 人
ちなみに，今月の支持率の 95% 信頼区間は，50.8% ～ 57.1%で，先月の支持率 50% はその範囲に含まれている。
先月のデータを探すのが面倒なので，回答者数が今回と同じ 978 人と仮定（標本は今月とは独立とする）
先月は，支持する:0.50*978 ≒ 489 人，支持しない:978-489 = 489 人
さて，「帰無仮説:今月と先月で内閣支持率は変化していない」を検定してみる。
> chisq.test(matrix(c(528, 450, 489, 489), 2))
Pearson's Chi-squared test with Yates' continuity correction
data: matrix(c(528, 450, 489, 489), 2)
X-squared = 2.9577, df = 1, p-value = 0.08547
内閣支持率には有意な変化があったとはいえない。
新聞のヘッドラインとはずいぶん違う。
新聞記事が，統計的検定結果を添えて記事を書いているのはかつて見たことがない。
また，マスコミ各社の内閣支持率は微妙に異なり，産経，読売などは高く，朝日，毎日などは低く出がちというのは，いったいどういうことだろう。
内閣支持者が，「済みません●●ですが，内閣支持率の調査をしています」という電話を受けたら，「がっちゃん!」と電話を切る。その逆もあり。
なんてことがあったのでは，正しい支持率が調査できるわけがない，ということではないだろうかなあ？？？
重要な世論調査については，実施主体者を明らかにしないという体制（というか公正な第三者機関が実施する）で行わなくてはならないのではないか？
実施主体で結果が変わるようでは統計学の名前が泣く。
かつてのアメリカの大統領選挙の予測と結果の乖離が教訓だろう。
追記:
2015年02月07日 読売新聞社は６～７日、全国世論調査を実施した。安倍内閣の支持率は５８％で、前回調査（１月９～１１日）の５３％から５ポイント上昇した。
調査の詳細はない（登録ユーザが見られるページは詳細が書かれているのかも知れない）
NHK の調査と同じ程度（978人）と仮定してみると，
> chisq.test(matrix(c(528, 450, 567, 411), 2))
Pearson's Chi-squared test with Yates' continuity correction
data: matrix(c(528, 450, 567, 411), 2)
X-squared = 2.9958, df = 1, p-value = 0.08348
NHK の結果と，読売新聞の結果は違っているとはいえないよ。
追記(2015/02/11)
「詳細情報記載ページ」を教えて戴いた。
電話番号を 5800 件用意し，うち有権者世帯が確認出来たのが 1913 件，世帯で有権者 1 人を無作為指定。有効回答は 1054 人（回答率 55%!! NHK の 65% より10% も低い!!「読売です」，「ガッチャン!!」なのか？？）
内閣支持率 58%（611 人），「支持以外」 443 人（しかし，「答えない」が 6% いるとしているのは，どんなものか？内閣支持率を求めるときは，分母から除くべきやいなや）
で，やり直し。
> chisq.test(matrix(c(528, 450, 611, 443), 2))
Pearson's Chi-squared test with Yates' continuity correction
data: matrix(c(528, 450, 611, 443), 2)
X-squared = 3.1056, df = 1, p-value = 0.07803
やはり，NHK の結果と違っているとは言えない。
==========
読売新聞が，うちの結果が正しいんだよというためには，調査対象数を倍増すればよい。
> chisq.test(matrix(c(528, 450, 567*2, 411*2), 2))
Pearson's Chi-squared test with Yates' continuity correction
data: matrix(c(528, 450, 567 * 2, 411 * 2), 2)
X-squared = 4.061, df = 1, p-value = 0.04388
59 文字解
f=function(s)s==paste(rev(strsplit(s,"")[[1]]),collapse="")
> f("abcあいう")
[1] FALSE
> f("abccba")
[1] TRUE
「時系列データと時間の相関について」だけど...
いわれているのはもっともだけど，そもそも対象とするようなデータの相関係数の分布が，普通の相関係数の分布とはまるっきり違うのだから，母相関係数の検定をそのまま使って結果がおかしいだろうという方が「おかしい」のでは？？
言いがかりというか，状況をわきまえない統計学の誤用。警鐘を鳴らす人も，そのあたりのことは分かっているんだろうかなあ。自分はわかっているけど，分かっていない人に警鐘を鳴らしているんだか，そんなことも分かっていなくて，おかしいよとだけいているんだか，よくわからない（自分はその人ではないのでね）。
sim = function(n=100) {
x = 1:n
y = cumsum(rnorm(n))
cor(x, y)
}
r = replicate(10000, sim())
hist(r)
つまりね。有意な相関があるという結果が出るのは「当たりまえに近い」ということ。
「いや，文字通り，当たり前に近いということをいっているのだよ！」ということなんだろうけどね。それこそ，「当たり前なんだよ」
ちょっとね，争点がフメイ。
------
念のため
一変量回帰分析（直線回帰）の傾きの検定と相関係数の検定は同じ...
treetreetreefiregoldtreetreetreetreegoldtreegoldtreegoldtreegoldtreegoldtreetreegoldtreegoldtreetreetreegoldtreetreetreewaterwatertreetreewatergoldgoldtreegoldtreegoldtreegoldtreegoldtreetreegoldtreegoldtreetreetreegoldtreetreefiregoldgoldwatermoontreemoontreegoldgoldtreewatergoldgoldtreetreegoldgoldtreetreemoonwatergoldgoldwatergoldtreetreemoongoldgoldmoon
を出力せよ
素直に書いて R で 171 文字解
策を弄して 166 154 150 文字解
GAWK では，BEGIN { } を除いて 178 165 文字解
当初の締め切り時間が過ぎたので，解をアップロードしておく
R は，以下の通り
a=gsub;cat(a(5,"moon",a(4,"water",a(3,"gold",a(2,"fire",a(1,"tree",a(6,31,"1112611166666166116114411436666616611612334515136436136154334615335")))))))
GAWK(MAWK）では，以下の通り（BEGIN{ と } を除く）
$0="1112611166666166116114411436666616611612334515136436136154334615335";gsub(6,31);gsub(1,"tree");gsub(2,"fire");gsub(3,"gold");gsub(4,"water");gsub(5,"moon");print
ちょっとしたことをやるだけなのに，わざわざパッケージを呼ぶなんて。
"How to Get the Frequency Table of a Categorical Variable as a Data Frame in R" だが，要するに，「度数分布表をデータフレームで返す関数 plyr:::count がいいよ！」ってことですな。
たとえば，iris データセットの Species の度数分布を知りたいとすると，
> ( w = table(iris$Species) )
setosa versicolor virginica
50 50 50
> ( t = data.frame(w) ) # as.data.frame でなくてよい
Var1 Freq
1 setosa 50
2 versicolor 50
3 virginica 50
> names(t)[1] = "Species" # 上の t の列名がイヤだというのだ
> t
Species Freq
1 setosa 50
2 versicolor 50
3 virginica 50
のようにしなきゃいけなくて，面倒！！ということのようだ。
plyr を使えば，簡単だよ！！！って。
> library(plyr)
> count(iris, 'Species') # なんで 列名を文字列で指定する？
Species freq
1 setosa 50
2 versicolor 50
3 virginica 50
おやおや。だなあ。
列名がイヤなら，table じゃなくて，xtabs を使えって，中澤さんも言ってるよ。
データフレーム名も情報に含んでくれるぞ（table 使いたきゃ，dnn を指定せよ）。
2 つくらいの関数を使うのを面倒というなら，関数化しておけ。
> data.frame(xtabs(~iris$Species))
iris.Species Freq
1 setosa 50
2 versicolor 50
3 virginica 50
「baloonplotの行名を並び替えたい」だが，細かいが気になることが満載だ
まず，日本語がだめ。「並び替え」ではなく「並べ替え」だ。
また，行名と列名について，間違って反対に覚えている。
次いで，R のプログラムだが，sample(c(1:16), 100, replace=TRUE) c(1:16) ってのが 1:16 と同じであることを分かっていない。それを直したとして，
sample(1:16, 100, replace=TRUE)は，遅い。以下のようにする。
sample(16, 100, replace=TRUE)とすべし。理由は sample のソースを見ればわかる。さすれば，
sample(c(0:25), 100, replace=TRUE)と同じ結果を得るための最適コードは，
sample(26, 100, replace=TRUE) - 1だな。
