Implemented type stable version of getproposal, fixes #12, confirmed …

…it results in faster algorithms

src/sampling.jl

@@ -1,14 +1,21 @@
 function getproposal(p::Prior, nparams)
 
  newparams = zeros(Float64, nparams)
-
- for i in 1:nparams
- newparams[i] = rand(p.distribution[i])
- end
+ update_newparams!(newparams, p)
 
  return newparams
 end
 
+update_newparams!(newparams,p::Prior) = update_newparams!(newparams, 1, p.distribution...)
+
+@inline function update_newparams!(newparams, i, x, y...)
+ newparams[i] = rand(x)
+ update_newparams!(newparams, i + 1, y...)
+end
+@inline function update_newparams!(newparams,i,x)
+ newparams[i] = rand(x)
+end
+
 particleperturbationkernel(x0, scale) = rand(Uniform(x0 - scale, x0 + scale))
 
 function perturbparticle(particle)
@@ -58,8 +65,6 @@ function getmodelfreq(particles, ABCsetup)
  return freq
 end
 
-
-
 function getmodelprob(currmodel, prevmodel, modelprob, ABCsetup)
 
  prob = ABCsetup.modelkern
@@ -74,7 +79,6 @@ function getmodelprob(currmodel, prevmodel, modelprob, ABCsetup)
 
 end
 
-
 function smcweights(particles, oldparticles, prior)
 
  weights = zeros(Float64, length(particles))

src/types.jl

@@ -6,8 +6,15 @@
 
  Create Prior type for ABC algorithm specifying priors for each parameters. This is an array of Distribution types from Distribution.jl, each elements corresponding to a parameter.
 """
+
+type Prior2
+ distribution
+ Prior2(distributionarray::Array{Distributions.Distribution{Distributions.Univariate,Distributions.Continuous},1}) = new(tuple(distributionarray...))
+end
+
 type Prior
- distribution::Array{Distributions.Distribution{Distributions.Univariate,Distributions.Continuous},1}
+ distribution
+ Prior(distributionarray) = new(tuple(distributionarray...))
 end
 
 type ParticleRejection <: Particle

test/test_bayesfactor.jl

@@ -65,9 +65,9 @@ nparticles = 500, maxiterations = 10^7)
 
 td, pM1 = generatedata()
 println("\t Checking ABC rejection model selection")
-abcresrej = runabc(ABCsetuprej, td);
+@time abcresrej = runabc(ABCsetuprej, td);
 @test isapprox(pM1, abcresrej.modelfreq[1], rtol = 0.05)
 
 println("\t Checking ABC SMC model selection")
-abcressmc = runabc(ABCsetupsmc, td);
+@time abcressmc = runabc(ABCsetupsmc, td);
 @test isapprox(pM1, abcressmc.modelprob[1], rtol = 0.05)

test/test_sampling.jl

@@ -62,3 +62,41 @@ weightsA, modelprob = ApproxBayes.getparticleweights(oldparticles, ABCsetup)
 for i in 1:ABCsetup.nmodels
  @test (map(x -> x.model, oldparticles).==i) == (weightsA[i, :].> 0.0)
 end
+
+
+
+#test get proposal function
+ABCsetup = ABCSMC(getnormal, 2, 0.1, Prior((Uniform(0.0, 20.0), Exponential(1.0))); nparticles = 100, maxiterations = 10^5)
+Niterations = 10^6
+p1 = zeros(Float64, Niterations)
+p2 = zeros(Float64, Niterations)
+@time for i in 1:Niterations
+ p1[i], p2[i] = ApproxBayes.getproposal(ABCsetup.prior, 2)
+end
+@test isapprox(mean(p1), 10.0, rtol = 0.001)
+@test isapprox(mean(p2), 1.0, rtol = 0.001)
+
+#test get proposal function with models with different Priors
+ABCsetup = ABCSMCModel([getnormal, getuniformdist], [2, 2], 1.0,
+[Prior([Exponential(1), Uniform(0.0, 1.0)]), Prior([Uniform(0.0, 3.0), Normal(2.0, 0.1)])]; nparticles = 100, maxiterations = 10^5)
+Niterations = 10^6
+p1 = Float64[]
+p2 = Float64[]
+p3 = Float64[]
+p4 = Float64[]
+@time for i in 1:Niterations
+ m = rand(1:2)
+ if m == 1
+ x = ApproxBayes.getproposal(ABCsetup.Models[m].prior, 2)
+ push!(p1, x[1])
+ push!(p2, x[2])
+ else
+ x = ApproxBayes.getproposal(ABCsetup.Models[m].prior, 2)
+ push!(p3, x[1])
+ push!(p4, x[2])
+ end
+end
+@test isapprox(mean(p1), 1.0, rtol = 0.01)
+@test isapprox(mean(p2), 0.5, rtol = 0.01)
+@test isapprox(mean(p3), 1.5, rtol = 0.01)
+@test isapprox(mean(p4), 2.0, rtol = 0.01)