For this kind of problem, you could allocate a big matrix in advance and fill
it in column-by-column, rather than push onto an array. Something roughly like
the following would work:
theta = Array(Float64, 3, 10_000)
theta[:, 1] = [2.0, 3.0, 4.0]
for k= 2:10_000
theta_cur = theta[:, k - 1]
theta_prop = theta_cur + rand(3)
likelihoodratio = computelr(theta_prop, theta_curr)
if rand() < min(1,likelihoodratio)
theta[:, k] = theta_prop
else
theta[:, k] = theta_curr
end
end
You could make something much faster still with devectorization. The code for
simulated annealing in Optim.jl might be helpful in this regard, although it
doesn't retain the full history of the chain.
-- John
On May 1, 2014, at 9:46 AM, Ethan Anderes <[email protected]> wrote:
> Ok, here is quickest example that I came up with.
>
>
> theta_init = [2.0, 3.0, 4.0]
> chain = Array{Float64,1}[init]
> for k=1:10_000
> theta_prop = chain[end]
> theta_prop += rand(3)
> likelihoodratio = computelr(theta_prop, theta_curr)
> if rand() < minimum([1,likelihoodratio])
> push!(chain, theta_prop)
> else
> push!(chain, theta_curr)
> end
> end
>
> I'm generating a markov chain of vectors. The code proposes a new state
> vector, theta_prop, and pushes it to chain if it satisfies a criterion. The
> problem is that I may run it for 10_000 steps, look at the results and then
> decide to run it for longer... so preinitializing isn't as natural.
>
> BTW: thanks a ton...I can't believe I'm having Stefan look at my code!
>
>
