裏 RjpWiki

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

変進小数の足し算

2016年08月10日 | ブログラミング
締め切りが 8/10 10.00 AM なのでその一分後に掲載されるように予約登録

変進小数の足し算
【概要】

小数の足し算をして下さい。ただしこの小数は、

小数第 n 位 が 11-n 進数

という不思議なルールになっています。例えば、小数第1位は 10進数、小数第9位は2進数です。
このルールを「変進小数」と呼びます。



【入出力】

入力は
8.622+3.177
こんな感じです。
2個の変進小数が「+」で区切られて並んでいます。

 

出力は、
11.811
のように、足し算の結果を変進小数で出力して下さい。


素直にプログラム

conv = function(s) {
    s = unlist(strsplit(s, "\\."))
    t = as.numeric(unlist(strsplit(s[2], "")))
    c(as.numeric(s[1]), c(t, rep(0, 9))[1:9])
}
deconv = function(x) {
    for (i in 10:2) {
        if (x[i] >= 12-i) {
            x[i-1] = x[i-1]+1
            x[i] = x[i]-(12-i)
        }
    }
    as.numeric(paste(x[1], ".", paste(x[2:10], collapse=""), sep=""))
}
x = "8.622+3.177"
options(scipen=100)
a = unlist(strsplit(x, "\\+"))
deconv(conv(a[1])+conv(a[2]))

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

念のため(^_^)

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

> 8月3日

> 10%水準だと強い主張ができないとお悩みのあなた!「片側検定」を使えば、10%水準の微妙な結果もあっという間に5%水準で統計的に有意になり、議論の説得性が増します!American Sociological Reviewでも使われている権威ある手法です!( ´ ▽ ` )ノ

皮肉っているのではあるが,この手法は片側検定がある場合しか使えない。

平均値の差の検定であれば,「2群の平均値の差の検定」などの場合だけに有効であるが,「3群以上の平均値の差の検定」(いわゆる一元配置分散分析)の場合には使えない。同じく,比率の差の検定も2群の場合のみ,独立性の検定ならば 2×2 分割表の場合のみなど。

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

信頼区間が重なっていても,有意差があることもあるんだよ

2016年07月28日 | 統計学

「有意差なしのロゴです。有意差がないときにご利用ください」だって...

https://pbs.twimg.com/media/CoTKh__VYAEhjkT.jpg

図が傾いているし,左右の棒の基線がなんなのか?

一番問題なのは,エラーバーがなんなのか。

たとえば,左側の棒の下端が第一群の標本平均,右側の棒の上端が第二群の標本平均として,エラーバーはそれぞれの母平均の95%信頼区間を表しているものとする。

以下のテストデータで,提示されロゴが表す状態が再現できる。

n1 = 40; mean1 = -39; se1 =   6/qt(0.975, n1-1); sd1 = sqrt(n1)*se1
n2 = 30; mean2 =  83; se2 = 128/qt(0.975, n2-1); sd2 = sqrt(n2)*se2

x = scale(rnorm(n1))*sd1+mean1
y = scale(rnorm(n2))*sd2+mean2
d = data.frame(g=rep(1:2, c(n1, n2)), z=c(x, y))

plot(c(0, 3), c(-50, 220), type="n", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
w = 0.15
rect(1-w, c(0, mean1), 1+w, c(mean1, 0))
rect(2-w, c(0, mean2), 2+w, c(mean2, 0))
abline(h=0, col="red", lty=3)
arrows(1:2, c(mean1-6, mean2-128), 1:2, c(mean1+6, mean2+128), angle=90, length=w*0.7, code=3)
text(1+1.2*w, mean1, paste("mean =", mean1), pos=4)
text(2+1.2*w, mean2, paste("mean =", mean2), pos=4)
text(2+1.2*w, 0, "base line", pos=4)
text(1:2+1.2*w, mean1-se1*qt(0.975, n1-1), "mean-se = -45", pos=4)

では,このデータで平均値の差の検定を行うとどうなるか

> t.test(x, y, var.equal=TRUE)

    Two Sample t-test

data:  x and y
t = -2.2519, df = 68, p-value = 0.02756

となり,信頼区間が重なっているにもかかわらず「有意な差である」という結論になる。

もっとも,var.equal=FALSE にすると,サンプルサイズをいろいろ変化させても有意な差ではないという結果にはなる。

 

> t.test(x, y, var.equal=FALSE)

    Welch Two Sample t-test

data:  x and y
t = -1.9472, df = 29.13, p-value = 0.06121

いずれにせよ,「エラーバーが重なっているから有意差はない」とはいい切れないのである。

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

ダミー変数をカテゴリー変数に変換

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

奥村先生が,名古屋市のHPVVのデータの解析をしている

kaito = read.csv("kaito.csv", header=FALSE, colClasses="character", fileEncoding="UTF-8")
dim(kaito)
birth = ifelse(rowSums(kaito[,5:11]=="1") != 1, NA,
  ifelse(kaito[,5]=="1", 6,
  ifelse(kaito[,6]=="1", 7,
  ifelse(kaito[,7]=="1", 8,
  ifelse(kaito[,8]=="1", 9,
  ifelse(kaito[,9]=="1", 10,
  ifelse(kaito[,10]=="1", 11,
  ifelse(kaito[,11]=="1", 12, NA))))))))
table(birth, useNA="ifany")

birth はご本人も「もっとエレガントな方法があるだろうが,面倒なので…」と書いてあるが,簡単なので以下のように

birth2 = as.matrix(kaito[,5:11] == 1) %*% 6:12
birth2 = ifelse(birth2 < 6 | birth2 > 12, NA, birth2)
table(birth2, useNA="ifany")

kaito[,5:11] == 1 は,元のデータを colClasses="character" で読んだので,必要悪
掛けるベクトルの選択には注意(6:12 なら問題なし)
数値で読んでいればそのまま行列掛算
無視すべき(NA にすべき)データは,ifelse(birth2 < 6 | birth2 > 12, NA, birth2) で処理

それにしても,ダミー変数展開をそのままデータ入力するとは非効率甚だしい

---------

質問:以下の項目で当てはまるものにチェックしてください

□ Aです □ Bです □ Cです □ Dです ✓ Eです □ Fです □ Gです

0, 0, 0, 1, 0, 0, 0 って入力したのか...バカ

普通は

質問:以下の項目で当てはまる番号に○をつけてください

1. Aです 2. Bです 3. Cです 4. Dです 5. Eです 6. Fです 7. Gです 8. いずれでもない 9. わからない

一桁の数字で入力するだろ

また,いずれでもない,わからない などの項目も設けていないと,無回答なのか,該当しないのか,単に回答し忘れたのかなどすら区別できない

統計学以前,統計調査票(アンケート調査票)の作り方を知らない素人がやった調査なのか?

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

アダムズ方式で議席数を計算

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

当初決められていた締め切りが 6月21日(火)10:00 AM なので,その 1 分後に投稿されるように予約
なお,締め切り日が延期されても,関知しない

衆院選挙制度改革法が成立し、「アダムズ方式」による議席数の割り当てが2020年の国勢調査に基づき用いられることが決まりました。
アダムズ方式は、各都道府県の人口を「ある同じ数値」で割って、その答えの合計が全国の議席数と同一になるように割る値を調整する計算方式です。
商が小数になる場合は切り上げることになっています。


例えば、250人、200人、150人の3つの県から10議席を選ぶとき、それぞれ65で割ると3.84…、3.07…、2.30…なので4,4,3となり合計が10になりません。
それぞれ75で割ると3.33…、2.66…、2なので4,3,2となりこれも合計が10になりません。
しかし、それぞれ70で割ると3.57…、2.85…、2.14…なので4,3,3となり合計が10になります。


標準入力の一行目に「議席数の合計」と「小選挙区の数」、二行目以降に小選挙区の数の分だけ人口が与えられます。
このとき、標準出力に各小選挙区の議席数を出力してください。


なお、人口として与えられる数は最大で1500万、小選挙区の数は最大で50とします。
「議席数の合計」「小選挙区の数」「人口」「割る数」のいずれも整数です。
※各小選挙区の議席数が一意に決まらないような入力は考えなくてもよいものとします。


【入出力サンプル】
標準入力
10,3
250
200
150

標準出力
4
3
3

例解

EPS = 1-1e-14

root = function(k, n) {
    f = function(x) sum(as.integer(k/x+EPS))-n
    left = 1
    right = max(k)
    repeat {
        mid = as.integer((left+right)/2)
        if (100 > right - left) {
            return(c(left, right))
        } else if (f(mid)*f(right) > 0) {
            right = mid
        } else {
            left = mid
        }
    }
}

func = function(s) {
    n = as.integer(unlist(strsplit(s[1], ",")))[1]
    k = as.numeric(s[-1])
    ans = root(k, n)
    for (i in ans[1]:ans[2]) {
        l = as.integer(k/i+EPS)
        if (sum(l) == n) {
            cat(paste(l, "\n", sep="", collapse=""))
            break
        }
    }
}

con = file("stdin", "r")
s = readLines(con)
func(s)

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

排他的 n 乗数(その2)

2016年06月15日 | ブログラミング

table は高くつくので以下のようにすれば十数倍速くなる。だが,まだまだ時間が掛かりすぎる。

func = function(s) {
    options(scipen=100)
    s = as.integer(unlist(strsplit(s, ",")))
    m = s[1]-1
    n = s[2]
    for (i in m:1) {
        a = as.integer(strsplit(as.character(i), "")[[1]])
        y = integer(10)
        y[a+1] = 1
        if (sum(y) == length(a)) {
            b = strsplit(as.character(i^n), "")[[1]]
            if (!any(b %in% a)) return(i)
        }
    }
    return("-")
}

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

排他的 n 乗数

2016年06月13日 | ブログラミング

「排他的 n 乗数」と呼ばれる数は,
    (a) 10 進数で表記したときに,同じ数字を含まない
    (b) n 乗した結果(10 進数表記)に,元の数に現れる数字が含まれない
というもので,例えば,608 は,608^3=224755712 なので「排他的3乗数」,284 は 284^4 = 6505390336 なので「排他的4乗数」である。

しかし,22^2=484 は,条件 (a) に合わないので 22 は「排他的 2 乗数」ではなく,123^4=228886641 は,条件 (b) に合わないので,123 は「排他的 4 乗数」ではない。

2 つの正の整数 m と n があるとき,m より小さい「排他的 n 乗数」のうち,もっとも大きな値を求めよ。そのようなものがない場合には "-" を返すものとする。

以下の R プログラムでは,時間が掛かりすぎる。

func = function(m, n) {
  options(scipen=100)
  for (i in (m-1):1) {
    a = unlist(strsplit(as.character(i), ""))
    if (!any(table(a) > 1)) {
      b = unlist(strsplit(as.character(i^n), ""))
      if (!any(b %in% a)) return(i)
    }
  }
  return("-")
}

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

スーパーマーズ?

2016年06月02日 | 雑感

> 31日、火星が地球に最接近します。「スーパーマーズ」とも言われ、

そんな言葉聞いたこともない

最接近とは言っても,当然,目で見えるほどの大きさになるわけでもなく,

数十パーセントも大きく見えるわけでもないので,「これぞ誇大広告!!」の最たるもんでしょう。

また,31日が過ぎると,突如普通の大きさに戻るわけでもないので,慌てる必要もない。

別のサイトで「どこで見られる?」ともあったが,「晴れていれば,日本国内,どこで見ても同じ」とうのが正解。

曇っていたら見えないでしょうが!!!!!

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

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

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重クロス集計をして確認しないと,因果関係は論じられない。

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