+ "content": "using ParetoSmooth\nusing Test\nusing Statistics\nusing AxisKeys\n\nimport RData\n\n\nlet og_array = RData.load(\"Example_Log_Likelihood_Array.RData\")[\"x\"]\n global log_lik_arr = copy(permutedims(og_array, [3, 1, 2]))\nend\nlet og_weights = RData.load(\"Weight_Matrix.RData\")[\"weightMatrix\"]\n global r_weights = exp.(permutedims(reshape(og_weights, 500, 2, 32), [3, 1, 2]))\nend\nr_eff = RData.load(\"Rel_Eff.RData\")[\"rel_eff\"]\nr_psis = RData.load(\"Psis_Object.RData\")[\"x\"]\nr_tail_len = Int.(RData.load(\"Tail_Vector.RData\")[\"tail\"])\npareto_k = RData.load(\"Pareto_K.RData\")[\"pareto_k\"]\nr_loo = RData.load(\"Example_Loo.RData\")[\"example_loo\"]\n\n\n# Add labels, reformat\nr_loo[\"pointwise\"] = KeyedArray(r_loo[\"pointwise\"][:, Not(4)];\n data = 1:size(r_loo[\"pointwise\"], 1),\n statistic=[:est_score, :mcse, :est_overfit, :pareto_k],\n )\n\nr_loo[\"estimates\"] = KeyedArray(r_loo[\"estimates\"];\n criterion=[:total_score, :overfit, :avg_score],\n estimate=[:Estimate, :SE],\n )\nr_loo[\"estimates\"](criterion=:avg_score) .= \n r_loo[\"estimates\"](criterion=:total_score) / size(r_loo[\"pointwise\"], 1)\n\n\n@testset \"ParetoSmooth.jl\" begin\n\n # All of these should run\n with_r_eff = psis(log_lik_arr, r_eff)\n jul_psis = psis(log_lik_arr)\n log_lik_mat = reshape(log_lik_arr, 32, 1000)\n chain_index = vcat(fill(1, 500), fill(2, 500))\n matrix_psis = psis(log_lik_mat; chain_index=chain_index)\n log_psis = psis(log_lik_arr; log_weights=true)\n\n jul_loo = loo(log_lik_arr)\n r_eff_loo = psis_loo(log_lik_arr, r_eff)\n\n display(jul_loo)\n \n # max 10% difference in tail length calc between Julia and R\n @test maximum(abs.(log.(jul_psis.tail_len ./ r_tail_len))) ≤ .1\n @test maximum(abs.(jul_psis.tail_len .- r_tail_len)) ≤ 10\n @test maximum(abs.(with_r_eff.tail_len .- r_tail_len)) ≤ 1\n \n # RMSE from R version is less than .1%\n @test sqrt(mean((with_r_eff.weights ./ r_weights .- 1).^2)) ≤ .001\n # RMSE less than .2% when using InferenceDiagnostics' ESS\n @test sqrt(mean((jul_psis.weights ./ r_weights .- 1).^2)) ≤ .002\n # Max difference is 1%\n @test maximum(log_psis.weights .- log.(r_weights)) ≤ .01\n\n\n ## Test difference in loo pointwise results\n\n # Different r_eff\n errs = (r_loo[\"pointwise\"] - jul_loo.pointwise).^2\n @test sqrt(mean(errs(:est_score))) ≤ .01\n @test sqrt(mean(errs(:est_overfit))) ≤ .01\n @test sqrt(mean(errs(:pareto_k))) ≤ .05\n errs_mcse = log.(r_loo[\"pointwise\"](:mcse) ./ jul_loo.pointwise(:mcse)).^2\n # @test sqrt(mean(errs_mcse)) ≤ .1\n\n # Same r_eff\n errs = (r_loo[\"pointwise\"] - r_eff_loo.pointwise).^2\n @test sqrt(mean(errs(:est_score))) ≤ .01\n @test sqrt(mean(errs(:est_overfit))) ≤ .01\n @test sqrt(mean(errs(:pareto_k))) ≤ .05\n errs_mcse = log.(r_loo[\"pointwise\"](:mcse) ./ r_eff_loo.pointwise(:mcse)).^2\n # @test sqrt(mean(errs_mcse)) ≤ .1\n \n # Test estimates\n errs = r_loo[\"estimates\"] - jul_loo.estimates\n @test maximum(abs.(errs)) ≤ .01\n \n errs = r_loo[\"estimates\"] - r_eff_loo.estimates\n @test maximum(abs.(errs)) ≤ .01\n\n # Test for calling correct method\n @test jul_loo.psis_object.weights ≈ psis(-log_lik_arr).weights\n @test r_eff_loo.psis_object.weights ≈ psis(-log_lik_arr, r_eff).weights\nend\n"
0 commit comments