Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 一休さんのとんち話にもあったが，「一日目は1粒，二日目は倍の2粒，三日目は更にその倍の4粒というように増やして行くとき，一日に 8000万粒以上になるのは何日目か」というもの。 プログラムしてはいけない。ただ，そのためには log10(2)≒0.3010 を記憶していなければならないが（大学受験生なら知っているだろう）。 80000000 ≦ 2^(n-1) 3+7/log10(2) ≦ n-1 27.3 ≦ n よって，28日目 なお，27日目までの累計も 80000000 を超える。 sum(2^0 + 2^1 + 2^2 + … + 2^(n-1)) = 2^n - 1

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 これは，物理法則なので，「グラフ用紙に点を打って，点の近くを通る直線を描いて，新規の x 値に対する y 値を読む」だけで答えが出る（「有効数字 1 桁で答えよ」なんだからなおさら）。 d = data.frame( x = sqrt(c(2e31, 2e32, 4e32, 5e32, 6e32, 8e32)), y = c(0.424e5, 1.34e5, 1.90e5, 2.12e5, 2.32e5, 2.68e5)) plot(y ~ x, d, xlim=c(5e15, 7e16), ylim=c(5e4, 60e4), xlab="sqrt(天体の質量[kg])", ylab="脱出速度[km/s]") ans = lm(y ~ x, d) abline(ans) new.x

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 （2）と似ているが str = "141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 言語によって，gsub の仕様が少し変わるが基本は同じ。 BEGIN { str = "141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 > 最初の2桁は 14 なので、アルファベット・インデックス表の 「o」 にあたり、カウントしません。 > 次の2桁は 15 なので、アルファベット・インデックス表の 「p」 にあたり、カウントします。 以下のプログラムは，要求されている仕様と異なる。以下のプログラムでは，14, 41, 15  のようにして二桁ずつ区切る。ははは。 この解は，βバージョン str = "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 outer は，実に使いでのある関数だ # 差が 6 の素数の組 n = 2:500 m = n[1-outer(n, n)] index = which(outer(m,  m, "-") == 6, arr.ind=TRUE) cat(paste(apply(matrix(m[index], ncol=2), 1, function(x) sprintf("(%i,%i)", x[2], x[1])), collapse=",")) # または cat(paste(mapply(function(x, y) sprintf("(%i,%i)", x, y), m[index[,2]], m[index[,1]]), collapse=",")) 出力 (5,11),(7,13),(

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 1:N までの数列から，2の倍数(2, 4, 6, ...)と3の倍数(3, 6, 9, ...)を取り除いた数列の和 func1 = function(N) { n = 1000000 x = 1:n x[1:(n%/%2)*2] = 0 x[1:(n%/%3)*3] = 0 x = x[x != 0] if (N > length(x)) return(NA) else return(sum(x[1:N])) } func2 = function(N) { m = N%/%2 s = 1 for (i in seq_len(m)) { s = s+12*i s = s-(i == m && N%%2==0)*(6*i-1) } s } func3 = function(N) { m

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 func2 <- function(n) { stopifnot(n <= 1388888888888885) # 桁数の同じ数を一群とする (5,6,7,8,9), (10,11,...,99), (100,101,...,999), (1000,1001,...,9999),... head <- 5 # 群の最初の数値 5, 10, 100, 1000,... kosuu <- 5 # 群に含まれる数値の個数 5, 90, 900, 9000,... keta <- 1  # 群に含まれる数値の桁数 1, 2, 3, 4,... end <- 0 # 群の最後の数値の1の位までの文字数（最初を1と数えて） repeat { begin <- end + 1 # 群の最初の数値の最

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 負の数は10の補数とするので，引き算もできる N <- 31 # 符号を含めた桁数 Add <- function(a, b) { carry <- function(s) { for (i in 1:N) { if (s[i] > 9) { s[i+1] <- s[i+1]+1 s[i] <- s[i]-10 } } return(s) } decomp <- function(s) { if (sign <- grepl("^-", s)) s <- sub("^-", "", s) s <- rev(tail(c(rep(0, 31), as.integer(unlist(strsplit(s, "")))), 31)) if (sign) { s <- 9-s s[1] <-

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 R だと invisible(sapply(1:10, function(n) cat(choose(n, 0:n), "\n"))) for を使う方が短いが for (i in 1:10) cat(choose(i, 0:i), "\n") choose なんていう関数を使うなということならば x <- 1; for (i in 1:10) cat(x <- c(x, 0)+c(0, x), "\n") いずれにしろ，Ruby よりは短いしトリッキーでもない

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 http://r-statistics-fan.hatenablog.com/entry/2014/08/09/103123 「Rで多倍長精度計算～とりあえず10000の階乗を計算してみる」であるが... 小数点以下も表示されてしまうので，整数部分だけを表示するように別解を。 先頭桁から区切っていくという仕様にした。別仕様は練習用。 library(Rmpfr) one = mpfr(1,200000) x = 10000 x = 100 fact.x = Reduce("*", c(one, 1:x)) print.mpfr = function(mpfr, n=100, m=10, Rev=FALSE) { num = sub("\\..+", "", formatMpfr(mpf

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 http://abrahamcow.hatenablog.com/entry/2014/09/11/024924 「時系列データの相関係数はあてにならない……のか？　教えて下さい」なんだけど... 私のコメントが気に触ることが多いようなのですが（特に悪意はないつもりなんですけど，すみませんね） 私は，経済学とか時系列についてはよく知らないのですが，「これは「見せかけの相関（擬似相関；spurious correlation）」の例だ」ということならば，偏相関係数を考えればよいのではないでしょうかね？？社会学などでは当たり前のように使われていると思うのですが。 > set.seed(1) > x = cumsum(rnorm(100)) > set.seed(2) > y = cumsu

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 以下の図は，1983年から1989年までの日経平均株価のグラフである。 「折れ線グラフは，変化を表すものだから，縦軸の原点は不要」というのは，いささか乱暴（目的によるのだろうが）ではないか。ちなみに，日経平均株価は間隔尺度では無く比尺度（ただし，0ということはあり得るので，純粋の比尺度では無いだろうが，株価0というのはあり得ない（あったらこわい）値なので）。 縦軸の目盛りがちゃんと描かれていれば，初期から終期までどれくらいの変化があったのかは，暗算すれば分かる（まあ，4倍くらいになったと言うことだが）。 だからといって，以下に示す2つのグラフが同じといえる訳ではないのは明らか。下の図は，「暗算しなくても目分量で4倍くらいかな」とわかる。逆に言えば，上の図は，うっかり見ていると（そんな奴

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 ごもっとも。しかし，その陳述の妥当性は「...」の内容によりますね。 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 5

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 新しい記事「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^

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 新しい記事「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

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 新しい記事「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

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 http://minato.sip21c.org/seek_pi_by_Archimedes.R 中澤さんの，「100ビット実数を使って円周率を10桁正しく出すには正何角形まで必要かを求めるためのRコード」について Rmpfr というパッケージの使い方がこれだけでわかり，ありがたかった。 もともと，計算時間がかかるというプログラムではなかったが，ちょっとした工夫で計算時間が短くなることがわかったので，記録しておこう。 末尾に追記あり。 中澤さんのオリジナルプログラム func1 = function() { library(Rmpfr) n = 6 p = sqrt(mpfr(3, 100)) / mpfr(2, 100) a = sqrt(mpfr(1, 100) - p ^ 2)

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 http://d.hatena.ne.jp/a_bicky/20140713/1405211663 R の named vector から subset を得る際の最悪計算量は O(n^2) の前半部分で，データフレームの factor 列の情報に基づき新たな列を算出するのに，names を使うと速いよと書いてあった。 > nrow(students) [1] 3000000 > class(students$grade) [1] "factor" > levels(students$grade) [1] "可"   "不可" "優"   "良" > scores 優   良   可 不可 30   20   10    0 という状況で，students$grade が，優，良，可，

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学 某所で 2, 5, 10, 19, 54, 224, 312, 616, 888, 977 より小さい数字の中で，素数が幾つあるかを求める Java プログラムを募集している。 Java で書けといわれて書けるかどうかはプログラマの技量を測るには適当か（適切ではない）と思うが，こんな計算量の小さな問題を解くなら，R の gmp を使えば朝飯前。 何を使おうと，最短時間で解答を得られれば，何の問題も無いと思うけど。 で，Java での解を募集している訳なので，これをこのまま Java に出来るわけも無い（うん？Java からも gmp ライブラリを参照できるかな？）ので，ネタバレにはならないだろう。 「より小さい数字の中で」とあるので，ちょっとまぬけな対応をした。 > library(g