add capability for dynamic models

Author: itsdfish

Examples/Delay_Discounting/Run_Delay_Discounting.jl

@@ -63,7 +63,7 @@ end
 #######################################################################################
 #                              Random Experiment
 #######################################################################################
-randomizer = Optimizer(;design_list, parm_list, data_list, model, approach=Randomize);
+randomizer = Optimizer(;design_list, parm_list, data_list, model, design_type=Randomize);
 design = randomizer.best_design
 new_data = [:random, 0, mean_post(randomizer)..., std_post(randomizer)...]
 push!(df, new_data)

Examples/Dynamic/Delay_Discounting.jl

@@ -0,0 +1,22 @@
+discount(t, κ) = 1/(1 + κ*t)
+
+function prob(κ, τ, t_ss, t_ll, r_ss, r_ll)
+    u_ll = r_ll * discount(t_ll, κ)
+    u_ss = r_ss * discount(t_ss, κ)
+    return 1/(1 + exp(-τ * (u_ll - u_ss)))
+end
+
+function loglike(κ, τ, t_ss, t_ll, r_ss, r_ll, data)
+    p = prob(κ, τ, t_ss, t_ll, r_ss, r_ll)
+    p = p == 1 ? 1 - eps() : p
+    p = p == 0 ? eps() : p
+    LL = data ? log(p) : log(1 - p)
+    # println(" choice ", data, " kappa ", κ, " tau ", τ, " t_ss ", t_ss,
+    #     " t_ll ", t_ll, " r_ss ", r_ss, " r_ll ", r_ll, " LL ", LL)
+    return LL
+end
+
+function simulate(κ, τ, t_ss, t_ll, r_ss, r_ll)
+    p = prob(κ, τ, t_ss, t_ll, r_ss, r_ll)
+    return rand() ≤ p ? true : false
+end

Examples/Dynamic/Run_Dynamic.jl

@@ -0,0 +1,89 @@
+#######################################################################################
+#                                 Load Packages
+#######################################################################################
+# set the working directory to the directory in which this file is contained
+cd(@__DIR__)
+# load the package manager
+using Pkg
+# activate the project environment
+Pkg.activate("../../")
+using Revise, AdaptiveDesignOptimization, Random, UtilityModels, Distributions
+include("Delay_Discounting.jl")
+#######################################################################################
+#                                  Define Model
+#######################################################################################
+Random.seed!(1204)
+prior = [Uniform(-5, 5), Uniform(-5, 50)]
+
+model = Model(;prior, loglike)
+
+parm_list = (κ = range(-5, 0, length=50) .|> x->10^x, 
+   τ = range(0, 5, length=11)[2:end])
+
+# parm_list = (κ = [.1], 
+#    τ = [.2,.5])
+
+design_list = (
+    t_ss = [0.0], 
+    t_ll = [0.43, 0.714, 1, 2, 3,
+        4.3, 6.44, 8.6, 10.8, 12.9,
+        17.2, 21.5, 26, 52, 104,
+        156, 260, 520], 
+    r_ss = 12.5:12.5:787.5,
+    r_ll = [800.0]
+)
+
+# design_list = (
+#     t_ss = [0.0], 
+#     t_ll = [5.0, 10.0], 
+#     r_ss = [12.0, 20.0],
+#     r_ll = [80.0]
+# )
+
+data_list = (choice=[true, false],)
+#######################################################################################
+#                              Simulate Experiment
+#######################################################################################
+using DataFrames
+true_parms = (κ=.12, τ=1.5)
+n_trials = 100
+optimizer = Optimizer(;design_list, parm_list, data_list, model);
+design = optimizer.best_design
+df = DataFrame(design=Symbol[], trial=Int[], mean_κ=Float64[], mean_τ=Float64[],
+    std_κ=Float64[], std_τ=Float64[])
+new_data = [:optimal, 0, mean_post(optimizer)..., std_post(optimizer)...]
+push!(df, new_data)
+
+for trial in 1:n_trials
+    data = simulate(true_parms..., design...)
+    design = update!(optimizer, data)
+    new_data = [:optimal, trial, mean_post(optimizer)..., std_post(optimizer)...]
+    push!(df, new_data)
+end
+#######################################################################################
+#                              Random Experiment
+#######################################################################################
+randomizer = Optimizer(;design_list, parm_list, data_list, model, approach=Randomize);
+design = randomizer.best_design
+new_data = [:random, 0, mean_post(randomizer)..., std_post(randomizer)...]
+push!(df, new_data)
+
+for trial in 1:n_trials
+    data = simulate(true_parms..., design...)
+    design = update!(randomizer, data)
+    new_data = [:random, trial, mean_post(randomizer)..., std_post(randomizer)...]
+    push!(df, new_data)
+end
+#######################################################################################
+#                                 Plot Results
+#######################################################################################
+using StatsPlots
+@df df plot(:trial, :mean_κ, xlabel="trial", ylabel="mean κ", group=:design, grid=false)
+hline!([true_parms.κ], label="true")
+
+@df df plot(:trial, :mean_τ, xlabel="trial", ylabel="mean τ", group=:design, grid=false)
+hline!([true_parms.τ], label="true")
+
+@df df plot(:trial, :std_κ, xlabel="trial", ylabel="σ of κ", grid=false, group=:design, ylims=(0,.3))
+
+@df df plot(:trial, :std_τ, xlabel="trial", ylabel="σ of τ", grid=false, group=:design, ylims=(0,2))

Examples/Monetary_Gambles/Run_Monetary_Gamble.jl

@@ -33,14 +33,13 @@ design_names = (:p1,:v1,:p2,:v2)
 design_list = (design_names,design_vals[1:100])
 
 data_list = (choice=[true, false],)
-
-optimizer = Optimizer(;design_list, parm_list, data_list, model)
 #######################################################################################
 #                              Simulate Experiment
 #######################################################################################
 using DataFrames
 true_parms = (δ=-1.0, β=1.0, γ=.7, θ=1.5)
 n_trials = 100
+optimizer = Optimizer(;design_list, parm_list, data_list, model)
 design = optimizer.best_design
 df = DataFrame(design=Symbol[], trial=Int[], mean_δ=Float64[], mean_β=Float64[],
     mean_γ=Float64[], mean_θ=Float64[], std_δ=Float64[], std_β=Float64[],
@@ -57,7 +56,7 @@ end
 #######################################################################################
 #                              Random Experiment
 #######################################################################################
-randomizer = Randomizer(;design_list, parm_list, data_list, model);
+randomizer = Optimizer(;design_list, parm_list, data_list, model, design_type=Randomize);
 design = randomizer.best_design
 new_data = [:random, 0, mean_post(randomizer)..., std_post(randomizer)...]
 push!(df, new_data)

src/AdaptiveDesignOptimization.jl

@@ -1,6 +1,7 @@
 module AdaptiveDesignOptimization
     using Distributions, Parameters, StatsFuns
     export Optimizer, Model, Optimize, Randomize
+    export Dynamic, Static
     export update!, mean_post, std_post, get_best_design
 
     include("structs.jl")

src/functions.jl

@@ -11,10 +11,14 @@ function prior_probs(prior::Nothing, parm_grid)
     return fill(1/length(parm_grid), size(parm_grid))
 end
 
-function loglikelihood(model::Model, design_grid, parm_grid, data_grid)
+function loglikelihood(model::Model, design_grid, parm_grid, data_grid, model_type, model_state)
     return loglikelihood(model.loglike, design_grid, parm_grid, data_grid)
 end
 
+function loglikelihood(model::Model, design_grid, parm_grid, data_grid, model_type::Dyn, model_state)
+    return loglikelihood(model.loglike, design_grid, parm_grid, data_grid, model_state)
+end
+
 function loglikelihood(loglike, design_grid, parm_grid, data_grid)
     LLs = zeros(length(parm_grid), length(design_grid), length(data_grid))
     for (d, data) in enumerate(data_grid)
@@ -27,6 +31,38 @@ function loglikelihood(loglike, design_grid, parm_grid, data_grid)
     return LLs
 end
 
+function loglikelihood(loglike, design_grid, parm_grid, data_grid, model_state)
+    LLs = zeros(length(parm_grid), length(design_grid), length(data_grid))
+    i = 0
+    for (d, data) in enumerate(data_grid)
+        for (k,design) in enumerate(design_grid)
+            for (p,parms) in enumerate(parm_grid)
+                i += 1
+                LLs[p,k,d] = loglike(parms..., design..., data..., model_state[i]) 
+            end
+        end
+    end
+    return LLs
+end
+
+function loglikelihood!(optimizer) 
+    @unpack model, log_like, design_grid, parm_grid, data_grid, model_state = optimizer
+    return loglikelihood!(model.loglike, log_like, design_grid, parm_grid, data_grid, model_state)
+end
+
+function loglikelihood!(loglike, log_like, design_grid, parm_grid, data_grid, model_state)
+    i = 0
+    for (d, data) in enumerate(data_grid)
+        for (k,design) in enumerate(design_grid)
+            for (p,parms) in enumerate(parm_grid)
+                i += 1
+                log_like[p,k,d] = loglike(parms..., design..., data..., model_state[i]) 
+            end
+        end
+    end
+    return nothing
+end
+
 function marginal_log_like!(optimizer)
     @unpack marg_log_like,log_like,log_post = optimizer
     marg_log_like .= marginal_log_like(log_post, log_like)
@@ -120,12 +156,36 @@ function update!(optimizer, data)
     return best_design
 end
 
+function update!(optimizer::Optimizer{A,MT}, data, args...; kwargs...) where {A,MT<:Dyn}
+    update_posterior!(optimizer, data)
+    update_states!(optimizer, data, args...; kwargs...)
+    loglikelihood!(optimizer)
+    marginal_log_like!(optimizer)
+    marginal_entropy!(optimizer)
+    conditional_entropy!(optimizer)
+    mutual_information!(optimizer)
+    best_design = find_best_design!(optimizer)
+    return best_design
+end
+
 function update!(optimizer::Optimizer{A}, data) where {A<:Rand}
     update_posterior!(optimizer, data)
     best_design = find_best_design!(optimizer)
     return best_design
 end
 
+function update_states!(optimizer, obs_data, args...; kwargs...)
+    @unpack model_state, data_grid, design_grid, parm_grid, update_state! = optimizer
+    for (d, data) in enumerate(data_grid)
+        for (k,design) in enumerate(design_grid)
+            for (p,parms) in enumerate(parm_grid)
+                update_state!(model_state, parms, design, data, obs_data, args...; kwargs...)
+            end
+        end
+    end
+    return nothing
+end
+
 function to_grid(vals::NamedTuple)
     k = keys(vals)
     v = product(vals...) |> collect
@@ -171,8 +231,13 @@ function std_post(optimizer)
     return std_post(post, parm_grid)
 end
 
-function create_state(T, dims, args...; kwargs...)
+function create_state(model_type::Dyn, T, dims, args...; kwargs...)
     state = fill(T(args...; kwargs...), dims)
-    state .= deepcopy.(state)
+    return state .= deepcopy.(state)
 end
 
+function create_state(model_type::Stat, T, dims, args...; kwargs...)
+    return nothing
+end
+
+