stat_test5.rb

require 'statsample'


def d(s)
        c = s.split('').select { |x| x == '1' }.size
	d = s.empty? ? 0 : c.to_f / s.length
        return (0.5 - d).abs, d
        
end

def seq2(n, c)
	l = []
	n1 = n
	c.times \
	{
		l << d(n.to_s(2))
		n = n.odd? ? n * 3 + 1 : n / 2
		break if (n < n1 || n == 1)
	}
	return l
end

def stat(l, t)
        return 0, 0, 0 if (l.empty? || t == 0) 
        t = t2 = 0
        l.each \
        {
		|x|
		t += x
		t2 += x ** 2
        }
        c = l.size
        a = t.to_f / c
        z = t2.to_f / c - a ** 2
        sd = Math.sqrt(z < 0 ? 0 : z)

	mx = l.max
	return a / t, sd / t, mx.to_f / t
end

def data(l)
	l.map \
	{
		|z|
		ns = z['n'].to_s(2)
		nl = ns.length
		
		d = seq2(z['n'], 3)
		
		asdm1 = stat(ns.split(/0+/).map { |x| x.length }, nl)
		
		
		[z['x'], nl, z['k'], z['cg'], z['cg'].to_f / nl, d, asdm1].flatten
	}
end

def mult(l, x, c, v)
	
	l2 = []
	l[0].size.times \
	{
		|i|
		t = c
		l.size.times \
		{
			|j|
			next if (x[j].nil?)
			t += l[j][i] * x[j]
		}
		l2 << [l[v][i], t]
	}
	return l2
	
end

l = File.open('db.txt').readlines
$mx = 20

l1 = []
l2 = []
8.times \
{
	|i|
	
	(20...100).each \
	{
		|j|
		z = Kernel.eval(l[i * 100 + j])
		next if (z['cg'] < $mx)

		z['x'] = j
		(j < 50) ? l1 << z : l2 << z
	}
}

$stderr.puts("l1# = #{l1.size}")
$stderr.puts("l2# = #{l2.size}")

l1 = data(l1).transpose
l2 = data(l2).sort_by { |x| x[1] }.transpose

v = ['x', 'nl', 'k', 'cg', 'h', 'd1', 'd2', 'd3', 'd4', 'd5', 'd6', 'a', 'sd', 'm']

h = {}
l1[v.index('k')].each { |x| h[x] = 0 if (h[x].nil?); h[x] += 1 }
$stderr.puts(h.inspect)

a = {}
v.each_with_index \
{
	|k, i|
	next if (['x', 'nl', 'k', 'cg'].member?(k))
	a[k] = l1[i].to_vector()
}

ds = a.to_dataset()

v1 = 'h'
r = Statsample::Regression.multiple(ds, v1)

$stderr.puts(r.summary)


x = []
r.coeffs.each \
{
	|k, z|
	x[v.index(k)] = z
}

l = mult(l2, x, r.constant, v.index(v1))
l.each_with_index \
{ 
	|x, i| 
	puts(([i, l2[v.index('x')][i], l2[v.index('nl')][i], l2[v.index('cg')][i], l2[v.index('nl')][i] * x[1]] + x).join("\t")) 
}