新聞までが騙されちょんWWW
オボッシーが立ち上げたページだという
ではあるが，3/31 21:49現在で，すでに
になっている。
オボッシー擁護派は喜んでいるようだが，早まりトチ狂ったエイプリルフールなんだろうね。
キーワード:小保方晴子，スタップ細胞，STAP，理研，捏造，不正論文
> 「そもそも、21とか14という少ない検体数では偶然の産物である場合が多すぎて、科学的にどうということは論じられない。そのため検定を行うことはほと んど意味をもたず、だから発表では検定の結果を示さなかったのかもしれない。今後の計画では150人に検体を増やしていくらしいが、少なくとも現段階での データは科学的意味を全く持たない」
そんなばかなことはない。
得られているサンプルサイズにおいて結果を判断するのが検定だ。
サンプルサイズが小さいということは，誤差も大きいということ。
その誤差よりも測定結果の差が大きいかどうかをみるのが検定だ。
「ミルクティーを作るとき，ミルクに紅茶をいれるのと，紅茶にミルクをいれるのがおなじかどうか」という話を知らないとはいわせないぜ。
http://analyze1842.blog133.fc2.com/blog-entry-36.html
まさに，「両側F検定（プププッ，変な名称を使うなよ）」だぜ。
> 150人のHLAを調べた結果が14人のときと全く同じであるとは考えられない。
サンプリングに偏りがないならば，150 人のときの該当数は 14 人のときの該当数の 150/14 倍になることが期待される。もっとも，その場合には，期待される該当数はおおむね 11 ずつ変化するので，誤差は大きいだろうけど。
もし，そうならないということなら，（最初の14人と後の何人かにおいて）サンプリングに偏りがあったということだ（どちらがズレているかはわからないが）。時間経過で傾向が変わったとか，複数の母集団を対象にした（背景因子が異なる集団からサンプリングされた）とか，理由はいろいろ。
さらに，もし該当の割合が大きく変わらないならば，サンプルサイズが大きくなるということは p 値は小さくなるということになる。
サンプルを2倍集めただけで，結論は大きく変わる。
結局の所は，ちゃんとした結果を出すために，どれくらいのサンプルを取らないといけないかを事前に評価してからデータを集める必要があるということ。
「データが増えるとどうなるかわからない」というのは，「適当にデータを集めて検定する」というのと，本質的に同じだ。
> fisher.test(matrix(c(1139, 1827, 24, 18), 2))
Fisher's Exact Test for Count Data
data: matrix(c(1139, 1827, 24, 18), 2)
p-value = 0.01625
> fisher.test(matrix(c(1139, 1827, 24*2, 18*2), 2))
Fisher's Exact Test for Count Data
data: matrix(c(1139, 1827, 24 * 2, 18 * 2), 2)
p-value = 0.000886
締め切りが 2016/03/24 10:00 AM なので，その 1 分後に公開されるように予約投稿
【概要】
年月日を指定します。曜日を計算してください。
ただし、地球ではなく、遠い昔、はるか彼方の銀河系にあった CodeIQ 星の「ラゲ暦」と呼ばれるカレンダーです。
日本で現在使われているカレンダー（グレゴリオ暦）との違いを下表にまとめました。
項目 ラゲ暦
曜日 t, u, v, w, x, y, z
一年の日数 345 日または344 日
一年の月数 11ヶ月 (A月～K月)
一ヶ月の日数
■平年の場合
A, B, D, F, H, I, K月は31日
C, E, G, J月は32日
■うるう年の場合
A, B, D, E, F, H, I, K月は31日
C, G, J月は32日
※ グレゴリオ暦と異なり、うるう年の方が一日短い
うるう年 ラゲ暦年が20で割り切れる年はうるう年
ただし、ラゲ暦年が80で割り切れる年は平年
ただし、ラゲ暦年が4000で割り切れる年はうるう年
なお、ラゲ暦の最初の日、1600年A月1日は、t曜日です。
【入力と出力】
入力は
1600.A.1
のような形で、標準入力から来ます。
ドット区切りで、左から順に 年月日 です。「年」は4～7桁、「日」は1～2桁です。
ラゲ暦は、1600年A月1日以降です。
CodeIQ星は栄華を誇り、このあと百万年ほど繁栄を続けます。
というわけで、入力 1600年A月1日から、1000000年K月31日 までの範囲です。
出力は、曜日です。
アルファベット小文字を一文字出力してください。
【入出力の例】
入力 出力
1600.A.1 t
1601.B.7 x
1602.C.14 v
1603.D.19 z
1604.E.25 w
1605.F.31 u
【補足】
不正な入力（"1700.A.11.B"、"Use your force, Luke"など）に対応する必要はありません。
グレゴリオ暦と同じく、月の最初の日は 1日 です。0日はありません。
171 文字以下にはできそうになかった
s=scan(file("stdin","r"),"",sep=".")
y=(v=as.integer)(s[1])
E=y%%20|!y%%80&y%%4000
cat(letters[(y*2-y%/%20+y%/%80-y%/%4000-4-E+cumsum(c(0,3,3,4,3,3+E,3,4,3,3,4))[utf8ToInt(s[2])-64]+v(s[3]))%%7+20])
締め切りが 2016/03/24 10:00 AM なので，その 1 分後に公開されるように予約投稿
【概要】
年月日を指定します。曜日を計算してください。
ただし、地球ではなく、遠い昔、はるか彼方の銀河系にあった CodeIQ 星の「キュラ暦」と呼ばれるカレンダーです。
項目 キュラ暦
曜日 t, u, v, w, x, y, z
一年の日数 345 日
一年の月数 11ヶ月 (A月～K月)
一ヶ月の日数 A, B, D, F, H, I, K月は31日
C, E, G, J月は32日
うるう年 ない
なお、キュラ暦500年の最初の日は、t曜日です。
【入力と出力】
入力は
500.A.1
のような形で、標準入力から来ます。
ドット区切りで、左から順に 年月日 です。「年」は3～4桁、「日」は1～2桁です。
キュラ暦は、500年A月1日に始まり、1599年J月26日で終わるので、入力はこの範囲の値になります。
出力は、曜日です。
アルファベット小文字を一文字出力してください。
【入出力の例】
入力 出力
500.A.1 t
501.B.4 u
502.E.7 v
503.J.17 w
504.K.28 z
【補足】
不正な入力（"123.Z.456"、"Force will be with you."など）に対応する必要はありません。
グレゴリオ暦と同じく、月の最初の日は 1日 です。0日はありません。
（R だと，？）71 文字以下にはできそうになかったので，短くする努力はしない。
f = function(s) {
ymd = unlist(strsplit(s, "\\."))
y = as.integer(ymd[1])
m = which(ymd[2] == LETTERS)
d = as.integer(ymd[3])
days = cumsum(c(0,31,31,32,31,32,31,32,31,31,32,31))
day = (y-500)*345+days[m]+d-1
letters[day%%7+20]
}
con = file("stdin", "r")
s = readLines(con)
cat(f(s))
締め切りが 3/21 10:00 AM なので，その 1 分後にアップされるように予約
下図のように数が並んでいます。
「cF」のようにマスを指定するので、それに対応する数（入力が「cF」なら111）を出力するプログラムを書いてください。
マスの位置を示す記号は、アルファベット2文字で、1文字目がa～zの小文字、2文字目がA～Zの大文字です。
出題者は，『「メモ化」するプログラムを書くのが普通？』とかいっているけど，漸化式があるので，for の方が簡単で時間も短い。
digits = 18
add = function(a, b, c) {
abc = a+b+c
carry = 0
for (i in seq_along(abc)) {
t = abc[i]+carry
abc[i] = t %% 10
carry = t %/% 10
}
abc
}
print.LongInt = function(a) {
if (all(a == 0)) {
a = 0
} else {
a = rev(a)
a = paste(rev(cutZero(rev(a[cumsum(a) > 0]))), collapse="")
}
cat(a, "\n", sep="")
}
n = 26
x = vector("list", n+1)
for (i in seq_along(x)) {
x[[i]] = vector("list", n+2)
}
zero = rep(0, digits)
one = zero
one[1] = 1
for (i in 1:(n+1)) for (j in 1:(n+2)) x[[i]][[j]] = zero
x[[1]][[3]] = one
for (i in 2:(n+1)) {
for (j in 3:(n+2)) {
x[[i]][[j]] = add(x[[i]][[j-2]], x[[i]][[j-1]], x[[i-1]][[j]])
}
}
#con = file("stdin", "r")
#line = readLines(con)
line = "zZ" # 188103751620198804
aB = unlist(strsplit(line, ""))
a = which(letters==aB[1])+1
B = which(LETTERS==aB[2])+2
print.LongInt(x[[a]][[B]])
--bignum を使える gawk なら，数列は以下の 8 行で定義できる
add3.awk:
BEGIN{
x[0, 1] = 1
for (i = 1; i <= 26; i++) {
for (j = 1; j <= 26; j++) {
x[i, j] = x[i, j-2] + x[i, j-1] + x[i-1, j]
}
}
}
$ gawk --bignum -f add3.awk
締め切りが 2016/03/15 10:00 AM なので，その 1 分後に投稿されるようにセット
クリスマスにプレゼントを配るサンタクロース。
配り始めると、すべて配るまで戻れません。
そこで、配るプレゼントをすべて持って出発します。
(帰ってくるときにはすべて配り終えているので、持っているプレゼントはなくなります。)
重たいプレゼントを持ったまま長距離を移動するのは大変ですので、サンタクロースも重たいプレゼントは先に配っていきたいと考えています。
そこで、できるだけサンタクロースの負担が少なくなるような順番を決めることにしました。
今回はこどもたちが用意した靴下の大きさが与えられます。
入れるプレゼントの重さは靴下の大きさと同じです。
（１の靴下には重さ１のプレゼント、2の靴下には重さ2のプレゼント、…）
サンタクロースはすべての靴下のサイズを知っており、靴下のサイズに合ったプレゼントを持って出発します。
訪問する家と家の間の距離は、それぞれの家にある靴下の大きさを掛けた値になっています。
例えば、サイズが2の靴下の家と、サイズが3の靴下の家の間の距離は6(2×3)です。
そして、この距離を移動するときには、「持っているプレゼントの重さ」と「距離」を掛けたエネルギーを消費します。
例えば、持っているプレゼントの重さが5で、家の間の距離が6のときは30(5×6)のエネルギーを消費します。
図のような靴下の家があった場合、
「1」→「2」→「3」→「4」の順にプレゼントを配ると、
9×2 ※9 ： 1を配ったので2,3,4のプレゼントの重さ 2 ： 1と2の距離
＋7×6 ※7 : 1と2を配ったので3,4のプレゼントの重さ 6 : 2と3の距離
＋4×12 ※1 : 1,2,3を配ったので4のプレゼントの重さ 12 : 3と4の距離
＝108
のエネルギーを消費します。
「4」→「3」→「2」→「1」の順にプレゼントを配ると、
6×12 ※6 ： 4を配ったので1,2,3のプレゼントの重さ 12 ： 4と3の距離
＋3×6 ※3 : 4と3を配ったので1,2のプレゼントの重さ 6 : 3と2の距離
＋1×2 ※1 : 4,3,2を配ったので1のプレゼントの重さ 2 : 2と1の距離
＝92
のエネルギーを消費します。
しかし、「4」→「1」→「3」→「2」の順にプレゼントを配ると、
6×4 ※6 : 4を配ったので1,2,3のプレゼントの重さ 4 : 4と1の距離
＋5×3 ※5 : 4,1を配ったので2,3のプレゼントの重さ 3 : 1と3の距離
＋2×6 ※2 : 4,3,1を配ったので2のプレゼントの重さ 6 : 3と2の距離
＝51
のエネルギーで済みます。
【問題】
標準入力から靴下の大きさが与えられます。（1行目にデータ数、2行目以降にそれぞれの家にある靴下の大きさ)
※靴下の大きさは正の整数で与えられ、すべて異なる大きさ、最大11個です。
同じ家には一度しか訪れることができず、必ずすべての家を訪問するとき、
必要なエネルギーが最小になる道順を求め、そのエネルギーを出力するプログラムを作成してください。
※図の場合は、この51というエネルギーが最小ですので、標準出力に51を出力します。
入力例
4
2
3
1
4
出力
51
=====
解答
最初に配るのは一番重いもの，
二番目に配るのは，一番軽いもの，
三番目に配るのは，残ったもので一番重いもの（最初の状態で二番目に重いもの）
四番目に配るのは，残ったもので一番軽いもの（最初の状態で二番目に軽いもの）
ということ（証明略(^_^)_）
ということで，以下のようなプログラムでいずれの例題も1秒以内で計算される。
f = function(x) {
x = sort(x)
n = length(x)
library(e1071)
if (5 >= n) {
p = permutations(n)
} else {
p = permutations(n-4)
y = x[-c(n, 1, n-1, 2)]
p = matrix(y[p], ncol=n-4)
p = cbind(x[n], x[1], x[n-1], x[2], p)
}
min.e = 1e200
junk = apply(p, 1, function(w) {
sum(sapply(2:n, function(i) sum(w[i:n])*w[i-1]*w[i]))
})
cat(min(junk))
}
system.time(f(c(2, 3, 1, 4))) # 51
system.time(f(1:11)) # 6112
system.time(f(c(20, 16, 14, 12, 9, 8, 5, 4, 3, 2, 1))) # 12004
> system.time(f(c(2, 3, 1, 4))) # 51
51 ユーザ システム 経過
0.016 0.001 0.011
> system.time(f(1:11)) # 6112
6112 ユーザ システム 経過
0.315 0.006 0.313
> system.time(f(c(20, 16, 14, 12, 9, 8, 5, 4, 3, 2, 1))) # 12004
12004 ユーザ システム 経過
0.300 0.003 0.295
R で，cor.test(x, y, method="spearman") とやれば，スピアマンの順位相関係数が表示されるが，信頼区間は表示されない。さて，どうしたものかというのがスタート地点。
上のように呼んだら，R は，x, y の代わりに rank(x), rank(y) を引数としてピアソンの積率相関係数を計算する。そして（理論的に）それは，スピアマンの順位相関係数に他ならない（実際，計算結果は一致する）。
では，cor.test(rank(x), rank(y)) としたらどうなるか。cor.test は x, y の測定値そのままではなくて，順位をデータだと思ってピアソンの積率相関係数を計算する。R は，計算されるのはスピアマンの順位相関係数だとわかっているからそのように計算する。しかし，ピアソンの積率相関係数の信頼区間は計算しない。
なぜ，計算しないか？順位をデータだと思ってピアソンの積率相関係数を計算するが，そのまま信頼区間の計算をしてよいかどうかについては「否」としているんだろう。
なぜか？元のデータは少なくとも正規分布するけど，順位に変換したデータは正規分布しないじゃないか。だったら，信頼区間は計算できないだろうということかな？
以下にシミュレーション・プログラムを示す。アイリス・データファイルからデータを取って，ブートストラップにより信頼区間を求めるプログラムである。（プログラムが間違っていたら元も子もないが）
当然ながら，両者が完全に一致するわけはないが，近似計算としては十分な意味を持つと思うがいかが？
CL.app というのは，データの代わりにデータの順位を与えてピアソンの積率相関係数の信頼区間を計算する方式で頸sんした信頼区間
CL.boot というのブートストラップによる信頼区間
まあ，結論は，どっちでもいいけどね（^_^;）
# ブートストラップ回数
N = 10000
# iris データセットの iris[1:50, 1:2] から，サンプルサイズ 30 で無作為抽出
j = sample(50, 30)
x = iris[j, 1]
y = iris[j, 2]
# rank(x), rank(y) で cor.test を呼べば，cor.test は順位相関係数を計算する
# おまけに，信頼区間も
a = cor.test(rank(x), rank(y))
a$estimate # スピアマンの順位相関係数（ただし，表示はピアソンの積率相関係数）
a$conf.int # 信頼区間
# ブートストラップ
z = sort(atanh(sapply(seq_len(N), function(i) {
i = sample(30, 30, replace=TRUE)
cor.test(x[i], y[i])$estimate
})))
tanh.z = tanh(z)
# ブートストラップによる信頼区間
boot = tanh.z[N*c(0.025, 0.975)]
layout(1:2)
par(mar=c(3, 3, 1.2, 0.5), mgp=c(1.8, 0.8, 0))
hist(z, breaks=30, probability=TRUE, main="Fisher's z transformation")
x = seq(0.2, 1.6, 0.01)
lines(x, dnorm(x, mean(z), sd(z)), col="red")
hist(tanh.z, breaks=30, probability=TRUE, main="Correlation coefficients")
legend("topleft", legend=c(
sprintf("Spearman's corr. =%.5f", a$estimate),
sprintf("CL.app =[%.5f, %.5f]", a$conf.int[1], a$conf.int[2]),
sprintf("CL.boot=[%.5f, %.5f]", boot[1], boot[2])),
bty="n")
layout(1)
現場では，こんな図を描いて教えているんだ。びっくり。
地球と月の大きさの割合と相互位置には正解を得られるように気を配ったのだろうが，観察値の地表から線を引くというのは，ばっかじゃないかと，皆が言う。
下のような，図を描くべき。
これだと月の高度が高すぎるけど，それは，月の位置をそこに指定してあるから。ちゃんとした高度が得られるように月の位置を変えないといけない。
金星の高度も高すぎるが，その原因は金星の公転軌道の半径が大きく描かれすぎているから。これも，ちゃんと直さないといけない。
この人，「『太陽に隠されるから見えない』なんておかしいだろ」っていってんだろうか。
東京都教育委員会もこの人もそうだけど，太陽が沈んだらすぐ星が見えるようになるって思っていない？
太陽と金星がこの図のような場合，真昼は当然としても，金星が太陽の東側にある場合は日没直後はまだ空が明るく，金星が太陽の西側にあるは場合は日の出前間近ではもう既に空が明るく，「肉眼では金星は見えない」。
図は（北極方向から見た図だとして），金星は太陽の東側にあるので，「明けの明星」の状態だけど，太陽から十分離れていないので，日の出前に明るくなってしまい，その状態では明けの明星としては観察できないだろう。
この人たちは，水星を見たことはないんだろうな。水星は，見かけでは太陽からあまり離れないから，なかなか見ることはできないのだよ。
なんか，頭の中の不完全な理論（？）で話をしているんじゃないか？
平成27年3月24日の東京都（139.73E, 35.70N）の日の入りは17:57のようで，「日没から30分後」は18:27ころ。
どんな夕空か？わかりますか？ 「天文薄明」という用語をご存じない？
まだ明るくて，月と金星以外は確認しづらい。火星はみえないかな。
火星より淡い，おひつじ座やうお座の星は，絶対に!!見えないだろう。
そもそも，その夕方は，晴れていたのかな？
.... 一応，晴れていたようだ。（って，そんな情報には何の意味もないのだけど）
http://weather.goo.ne.jp/past/days/03/24/662/
インターネットの世の中は，怖い....
なお，揚げ足取りに等しいが，夕空の惑星と月の配置関係も出題文とはだいぶ違う。
「日時を特定して出題したのは，リアリティーを出すため」と言っているが，それは全く逆。
出題する立場とすれば，特定の日時を設定して，その時にどのように見えて，実際の宇宙空間での位置関係を確かめた上で選択肢を作るべき，
本来は，出題者が実際に観察してそれをもとに出題すべき。
そして「問題文には，特定の日時など記述すべきではない」。
だって，特定の日時が記述されていても，解答時にインターネットを利用できない受験生には何のメリットもないし，不要の情報（リアリティーなんか，そんな場合には，必要ない）。
もし特定の日時を指定していていればその時の実際の惑星配置は特定できるわけで，それが間違えていれば明らかに誤った出題だから。言い逃れできない。
当然のことだけど，そのような状況ではこのような観察結果が得られるという論理的な思考過程があるべきだということはいうまでもない。
当日の金星の形は，問題文の図2とは全く異なる。
太陽と火星と金星の位置関係は
誰が調べても，こうなるんだから，正解例が不適切なのはあきらか。
言い逃れできない。
東京都立高校の理科の入試問題が物議を醸している
正解を「イ」としている
「ウ」じゃないの？？？ということで，「わーわーわーわー」（^_^;）いってます
締め切りが 2016/03/03 10:00 AM なので，10:01 に登録されるよう予約
自然数 n に対して、関数 F(n) を、nn（n の n 乗）を 10 進数で表したときの桁数と定義します。
例えば、55 = 3125 ですので、F(5) = 4 です。同様に、F(20) = 27 です。
さらに、2 以上の自然数 m に対して、F(n) = m となる n の値を G(m) と定義します。
もしそのような n が存在しない場合は、G(m) = -1 と定義します。
例えば G(4) = 5、G(27) = 20 です。同様に、G(7) = -1、G(101) = 57、G(11111) = 3173 となることが確かめられます。
標準入力から、自然数 m（2 ≦ m ≦ 1010）が与えられます。
標準出力に G(m) の値を出力するプログラムを書いてください。
チョー簡単
G = function(m) {
a = floor(uniroot(function(x) x*log10(x)-m,c(1, 1e10))$root)
if (ceiling(a*log10(a)) == m) cat(a) else cat(-1)
}
答が合わないときがあるんじゃないかというコメントをいただきました。
n が 10 の羃乗である場合は，間違った答になりますね。10^10 は 11 桁ですからね。
G = function(m) {
a = floor(uniroot(function(x) x*log10(x)-m,c(1, 1e10))$root)
b = a*log10(a)
if (b == m-1) return(a) else if (b == m || ceiling(a*log10(a)) != m) return(-1) else return(a)
}
でよいでしょうか。
締め切りの件は，先方が勝手に繰り延べるだけなので，こちらでは対応しません。
締め切りは 03/03 10:00AM なので，その 1 分後にこの回答が投稿されるように設定しておく。まだ 3 ヶ月以上も先のことで，鬼が笑うどころか，鬼が飽きてあくびをしたり，居眠りしそうなのだけど，しょうがない。
累乗の累乗，つまり，x^y^z = x^(y^z) は，x, y, z (>1) がちょっと大きくなると，とてつもなく大きな値になり，通常のコンピュータの演算ではオーバーフローして Inf という表示がされることになる。
x = 4.5
y = 6.7
z = 9.8
のとき，
> x^y^z
[1] Inf
x = 2.1
y = 7.2
z = 3.6
のときも，
> x^y^z
[1] Inf
となり，どっちが大きいかわからないのですけど？ということか。
この問題は，何組かの x, y, z が与えられたとき x^y^z が最も大きいのはどれかを判定せよというもの。
a = x^y^z の両辺の対数をとって
log(a) = y^z * log(x)
更にもう一回両辺の対数をとって
log(log(a)) = z * log(y) + log(log(x))
一般に，u < v なら log(u) < log(v) なので，
z * log(y) + log(log(x)) が大きい方が log(log(z * log(y) + log(log(x)))) も大きい。
log(log(z * log(y) + log(log(x)))) はコンピュータで表せないほど大きい値でも，z * log(y) + log(log(x)) はそんなに大きい値ではないことがある。
ということで，回答は log(log(z * log(y) + log(log(x)))) が最も大きいものが x^y^z も最も大きいということである（高校生レベルで解ける問題だ）。
コンピュータで表せる範囲内だと x^y^z は，exp(exp( z * log(y) + log(log(x)) )) であることは以下のように確かめることができる。
x = 1.2
y = 2.3
z = 3.4
> x^y^z
[1] 22.09543
> exp(exp( z * log(y) + log(log(x)) ))
[1] 22.09543
こんな簡単な問題を出す人もいるのだから，それに答える人がいてもよいではないかwww
追記 2021/06/05
Julia では
x1 = big"4.5" # 4.5
y1 = big"6.7" # 6.699999999999999999999999999999999999999999999999999999999999999999999999999972
z1 = big"9.8" # 9.800000000000000000000000000000000000000000000000000000000000000000000000000028
x1^y1^z1 # 8.141852291891192175294360390819933153170222320673963475352901403217459199966931e+81393094
x2 = big"2.1" # 2.099999999999999999999999999999999999999999999999999999999999999999999999999986
y2 = big"7.2" # 7.199999999999999999999999999999999999999999999999999999999999999999999999999972
z2 = big"3.6" # 3.599999999999999999999999999999999999999999999999999999999999999999999999999986
x2^y2^z2 # 1.384119872013753626430173206328153257960203418740675402484696742668262448719525e+393
z1 * log(y1) + log(log(x1)) # 19.0488334435159380073817381785428805350652170271541018037443913083315619703749
z2 * log(y2) + log(log(x2)) # 6.8082012132336013423449913855664088058111210532323767709225454899520398965019
log(log(z1 * log(y1) + log(log(x1)))) # 1.080789693280837499208500181126913735357567639885908188155661700356977982900217
log(log(z2 * log(y2) + log(log(x2)))) # 0.6513496823931141280708339598949722274473422067642296361862125798228163600446146
