johnmyleswhite

total = 0
N = 300

start_time = time()

for a in 0:(N - 1)
	for b in 0:(N - 1)
		for c in 0:(N - 1)
			if a^2 + b^2 == c^2
				total = total + 1

johnmyleswhite

using Distributions

d = Dirichlet([100.0, 17.0, 31.0, 45.0])

X = rand(d, 1_000_000)

fixed_point(X)
@elapsed alpha = fixed_point(X)
norm(d.alpha - alpha, Inf)

johnmyleswhite

# Implementation of fixed point algorithm described in
#  "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000
function invdigamma(y::Float64)
    # Closed form initial estimates
    if y >= -2.22
        x_old = exp(y) + 0.5
        x_new = x_old
    else
        x_old = -1.0 / (y - digamma(1.0))
        x_new = x_old

johnmyleswhite

fit(TDLearner, choices)
fit(HMM, signal)
fit(Lasso, Y ~ X, data, lambda = 1.0)

loglikelihood(model)

AIC(model)
BIC(model)

johnmyleswhite

# Uniform

a <- runif(10000, 0, 1)
b <- runif(10000, 0, 1)
c <- runif(10000, 0, 1)

df <- cbind(a, b, c)

df <- transform(df, s = a + b + c)

johnmyleswhite

function repeat{T}(v::Array{T},
                   dims::Integer...;
                   repetitions::Integer = 1)
    n = length(v)
    n_res = n * repetitions * prod(dims)
    res = Array(T, tuple(repetitions * n * dims[1], dims[2:end]...)...)
    for i in 0:(n_res - 1)
        index = mod(fld(i, repetitions), n)
        res[i + 1] = v[index + 1]
    end

johnmyleswhite

function bfs(start_node::Any,
			 next_nodes::Function,
	         at_node::Function)
	to_process = Array(Any, 0)
	depths = Array(Int, 0)
 
	push!(to_process, start_node)
	push!(depths, 0)
 
	while !isempty(to_process)

johnmyleswhite

# K clusters in N-dimensional space
# TODO: Fill out Distribution methods
# TODO: Restore empty cluster pruning
# TODO: Add method for getting point assignments (for all mixtures?)
# TODO: Make custom cont. univ. mixture type?
immutable MixtureMultivariateNormals <: ContinuousMultivariateDistribution
	mu::Matrix{Float64}
	sigma::Array{Float64}
	p::Vector{Float64}
    normals::Vector{MultivariateNormal}

johnmyleswhite

function betainc(x, pin, qin)
    # Are these appropriate?
    eps1 = eps(0.5)
    alneps = log(eps1)
    sml = eps(0.0)
    alnsml = log(sml)

    if x < 0.0 || x > 1.0
        error("x must lie in (0, 1)")
    end

johnmyleswhite

type PrimeArray
	n::BigInt
end

PrimeArray() = PrimeArray(BigInt(1))

function Base.getindex(a::PrimeArray, i::Integer)
	if i < 0 || i > length(Base.PRIMES)
		throw(BoundsError())
	end