裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

基本的な機能は,「基本的な重要性」を持っているのです

2016年05月30日 | ブログラミング

> みんなのRはテキストとしてかなり満足しているけど、str()・source()・write.table()・transform()あたりが載っていなくてびびった。

列挙に「・」を使う意味がわからないけど(^_^)

str() は,既存の関数の戻り値の内容ならば,"? 関数" を読めばわかるわけ

source() は,自作関数を作って活用する段階で必要なもの。それ以前なら,作り置いたファイルから,コンソールにコピペすれば十分なわけ

write.table() も,必要ではあるものの,read.table() したデータテーブルをいろいろ試行錯誤的というか追加的というか発展させていくわけで,write.table()したものを,そのまま使って分析へというのも,あんまり考えられない。データ入力からあらたな変数生成という過程をプログラムに残すということならば,むしろ write.table() で,ある時点のデータフレームを固定化する必要はないのではないか?むしろ,柔軟性を残す意味で,read.table() のあとの様々なデータ処理を一体化して保存しておくのがよいのではないか(このような時点では source() を使う必要もないのだよ)

transform() を便利だと思う人もいるだろうけど,直接的な記述をした方が「間違いがない」というメリットもあるんだよね。「わかりやすい」というのは,最大の長所。うっかりミスも防げる。

何が不要/必要,便利/不便というのは,それを使う人の技術水準で異なってくるんだね

指導者はある程度高度の水準にいるので,被指導者のレベルで何が必要かがわからなくなっているのではないか????

コメント
この記事をはてなブックマークに追加

正しい方法を使おう

2016年05月14日 | ブログラミング

2016 年 05 月 13 日 (金)

> R で計算するなら 1 - pnorm(60, 50, 10) あるいは pnorm(60, 50, 10, lower.tail = FALSE) twitter.com/yuba/status/73…

q が通常範囲であればどちらでもよいが,q が大きい場合は違いが生じる。

> options(digits=16)
> q = 100; 1 - pnorm(q, 50, 10); pnorm(q, 50, 10, lower.tail = FALSE)
[1] 2.866515719235352e-07
[1] 2.866515718791939e-07
> q = 120; 1 - pnorm(q, 50, 10); pnorm(q, 50, 10, lower.tail = FALSE)
[1] 1.27986510278788e-12
[1] 1.279812543885835e-12
> q = 150; 1 - pnorm(q, 50, 10); pnorm(q, 50, 10, lower.tail = FALSE)
[1] 0
[1] 7.619853024160527e-24


「統計学: R を用いた入門書 改訂第 2 版」にも,lower.tail を指定しない方法が書かれているが,どちらが正しいかは 3 番目の結果を見れば自明。

まあ,実際の統計学においては,その違いは結果を左右するようなものではないが。

コメント
この記事をはてなブックマークに追加

ほとんどの問題の締め切りが12月31日である

2016年05月10日 | 雑感

出題されている問題のうち,ほとんどすべての締め切りが12月31日である。

異常でしょう。

そもそも,このサイトの異常性(締め切りは常に延期される)からいえば,「何年の」12月31日であるかさえわからないという(嗤うしかない。「笑う」じゃないよ)。まあ,土曜日というのだから今年のことなんだろうけど。

実際の解答例や,出題者側からの例解発表もないしねえぇ。やってても,達成感がないのよねぇ。回答してみようという気が失せる。

コメント
この記事をはてなブックマークに追加

文政 7 年の算額

2016年05月05日 | ブログラミング

その1 外円(黒)の径が 1 尺のとき内円(赤)の径はいくらか(青は正六角形)

その2 大円(黒)の径が 3 尺 6 寸,中円(青)の径が 9 寸のとき,小円(赤)の径はいくらか

その3 直角三角形の面積が 756 歩,絃寸と股寸の和が 1 丈 4 尺 7 寸のとき,釣,股,絃はいくらか

 

注1 その1,その2 における「径」は「半径」

注2 1 丈 = 10 尺 = 100 寸

注3 その3における「歩」は 一辺が 1 寸の正方形の面積

コメント
この記事をはてなブックマークに追加

そんなもの,いらない!

2016年04月05日 | ブログラミング

> Surface Pro 4とSurface Book、岡山県備前市ふるさと納税の返礼

そんなもの,いらない! めいわく。

コメント
この記事をはてなブックマークに追加

ふーん。どういう肩書だ?

2016年04月04日 | 雑感

株式会社***** CEO, ****-****/**** ****店長, ****大学特任講師, **大学/***大学非常勤講師, はかせ(工学)

コメント
この記事をはてなブックマークに追加

炎上参加者の特性

2016年04月03日 | 統計学

http://mainichi.jp/articles/20160402/mog/00m/300/003000c

> (1)子供と同居している親は、そうでない人よりも、炎上行為に参加しやすい(2)個人年収や世帯年収が高くなるほど、炎上に参加する確率が高まる−−という結果が示された。「子供を持つ、裕福な人ほど炎上に参加している」(山口氏)ということになり

この研究,詳しくは書かれていないが,炎上に参加したかしないかを説明変数でロジスティック回帰したものと思われる。

(1) と (2) は,偏回帰係数を読み取っての解釈であろうが,そもそも偏回帰係数は他の説明変数が一定の場合の係数(だからこそ「偏」)なのだから,両者を併せて「「子供を持つ、裕福な人ほど炎上に参加」というのはおかしい。子供の有無と裕福かどうか,炎上参加・不参加の3重クロス集計をして確認しないと,因果関係は論じられない。

コメント
この記事をはてなブックマークに追加

一日早いエープリルフールか(その2)

2016年04月02日 | 雑感

https://stap-hope-page.com

404 Error になっているという情報もあったけど,04/02 18:30 現在,復活しているなぁ

彼女自身,このレシピで STAP 細胞を再現できなかったんじゃなかったっけ?

どうしようというんだろうか。

コメント
この記事をはてなブックマークに追加

一日早いエープリルフールか

2016年03月31日 | 雑感

新聞までが騙されちょんWWW

オボッシーが立ち上げたページだという

https://stap-hope-page.com

ではあるが,3/31 21:49現在で,すでに

403 - Forbidden Error

You are not allowed to access this address.
If the error persists, please contact the website webmaster.

になっている。

オボッシー擁護派は喜んでいるようだが,早まりトチ狂ったエイプリルフールなんだろうね。

キーワード:小保方晴子,スタップ細胞,STAP,理研,捏造,不正論文

コメント
この記事をはてなブックマークに追加

例数が少ないと検定はできないのか?

2016年03月24日 | 統計学

> 「そもそも、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日 | ブログラミング

締め切りが 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日 | ブログラミング

締め切りが 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))

コメント
この記事をはてなブックマークに追加

二次元の数列

2016年03月21日 | ブログラミング

締め切りが 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日 | ブログラミング

締め切りが 2016/03/15 10:00 AM なので,その 1 分後に投稿されるようにセット

クリスマスにプレゼントを配るサンタクロース。
配り始めると、すべて配るまで戻れません。
そこで、配るプレゼントをすべて持って出発します。
(帰ってくるときにはすべて配り終えているので、持っているプレゼントはなくなります。)

重たいプレゼントを持ったまま長距離を移動するのは大変ですので、サンタクロースも重たいプレゼントは先に配っていきたいと考えています。
そこで、できるだけサンタクロースの負担が少なくなるような順番を決めることにしました。

今回はこどもたちが用意した靴下の大きさが与えられます。
入れるプレゼントの重さは靴下の大きさと同じです。
(1の靴下には重さ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

コメント
この記事をはてなブックマークに追加

スピアマンの順位相関係数の信頼区間

2016年03月14日 | ブログラミング

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)

コメント
この記事をはてなブックマークに追加