Julia言語のお勉強をしているのだけれど, 多倍長精度浮動小数点数はBigFloatで直ぐに使えるし, Uniform Sum DistributionのMonte Carloは直ぐに出来てNapier数と戯れられるし, 頻度主義とベイズ法のブートストラップも容易に切り替えられる.
頻度主義:
using Distributions, OnlineStats
function bootCI(data, stat::Function, CI::Float64, reps::Integer)
tmp = similar(data)
boot = [stat(rand!(tmp, data)) for i in 1:reps]
low, high = quantile(boot, [(1-CI)/2, (1+CI)/2])
(value=stat(data), low=low, high=high)
end
ベイズ法:
using StatsBase, Distributions
function bayesbootCI(data, stat, CI::Float64, reps::Integer)
d = Dirichlet(length(data), 1)
boot = [stat(data, weights(rand(d))) for i in 1:reps]
low, high = quantile(boot, [(1-CI)/2, (1+CI)/2])
(value=stat(data), low=low, high=high)
end
Schemeであれだけ苦労したフラクタル図形の描写も, Juliaだとあっさり出来る. これはJulia集合.
さらには数独クイズの解法もあっさりプログラム出来る. そろそろ交換部品が無くなっているのじゃないかと思われる私の6年以上前の古いラップトップでも9☓9が0.05-0.5 secくらいですんなり解ける.
blockvalid(x, v) = count(isequal(v), x) ≤ 1
function backtrack!(x, z, idx)
idx > length(z) && return true
pos = z[idx]
iloc = 3div(pos[1]-1, 3)
jloc = 3div(pos[2]-1, 3)
filled = 0
@inbounds for k in 1:9
filled |= x[pos[1], k] | x[k, pos[2]]
end
@inbounds for k1 in 1:3, k2 in 1:3
filled |= x[iloc+k1, jloc+k2]
end
@inbounds for i in 1:9
k = 1 <<i
if k & filled == 0
x[pos] = k
backtrack!(x, z, idx+1) && return true
end
end
x[pos] = 1
return false
end
function ssolve(lines, i)
t = [1 << (lines[10i-j][k] - '0') for j in 8:-1:0, k in 1:9]
z = findall(isequal(1), t)
backtrack!(t, z, 1)
sum([100, 10, 1] .* trailing_zeros.(t[1, 1:3]))
end
lines = readlines("p096_sudoku.txt")
@time sum(ssolve(lines, i) for i in 1:50)
@time sum(ssolve(lines, i) for i in 1:50)
チンパンジーに画像への反応速度を応用したゲームを学習させ, ヒトと対戦させてもヒトは最早チンプの圧倒的速度に勝てないのと同様に, 計算機はちゃんとした方法を覚えこませたら圧倒的なまでの能力を発揮するのを実感した. モデルのプログラムは計算速度にあまり影響しなくて既存の統計的手法を使う所はRCallに任せて, 基本はJuliaを使おうと思う. その前にRやBioconductorで何が出来てどこから先は厳しいのかをきちんと調べ, またRStanでのベイズ統計モデリングをもう少し学習してから本業に取り組みたい.
【追記】Juliaで低水準言語による表現をチェックすれば訳の分からないことをしているのは一目瞭然だし, 変更可能な「構造体」mutableは数学的に使いやすそう.
