https://github.com/nwbarendregt/AdaptNormThresh
Raw File
Tip revision: 2878a3d9f5a3b9b89a0084a897bef3414e9de4a2 authored by nwbarendregt on 03 May 2022, 15:08:29 UTC
Add files via upload
Tip revision: 2878a3d
SC_Bellmans.m
% SC_Bellmans.m
% Function used to obtain normative thresholds via dynamic programming for 
% SNR change task from Barendregt et al., 2022.

function [theta,rho] = SC_Bellmans(T,dt,t_i,dg,m,c,R)

% Define time and likelihood discretizations:
t = 0:dt:T;
g = dg:dg:(1-dg); g_i = find(g==0.5);

% Calculate Gaussian means from SNR timeseries (assuming standard 
% deviation of 1):
mu = sqrt(unique(m)/2);

% Construct likelihood transfer functions given by Eq. (14):
P_gg = NaN(length(g),length(g),length(mu));
for i = 1:length(mu)
    P_gg(:,:,i) = diag(g)*exp(-0.5/dt*(log(((g-1)'*g)./(g'*(g-1)))/(2*mu(i))-mu(i)*dt).^2)/sqrt(2*pi*dt)./(2*g*mu(i)-2*g.^2*mu(i))+...
        diag(1-g)*exp(-0.5/dt*(log(((g-1)'*g)./(g'*(g-1)))/(2*mu(i))+mu(i)*dt).^2)/sqrt(2*pi*dt)./(2*g*mu(i)-2*g.^2*mu(i));
    P_gg(:,:,i) = P_gg(:,:,i)*diag(1./sum(P_gg(:,:,i),2));
end

% Initialize secant method:
rho = 0.1; tol = 1e-5; k = 1;

% Pre-allocate value function V and maximal index of value function V_I:
V = NaN(length(g),length(t)); V_I = NaN(length(g),length(t));

% Calculate value function using backward induction:
[V(:,end),V_I(:,end)] = max([g'*R-t_i*rho(k) (1-g)'*R-t_i*rho(k)],[],2);
for j = (length(t)-1):-1:1
    [V(:,j),V_I(:,j)] = max([g'*R-t_i*rho(k) (1-g)'*R-t_i*rho(k) squeeze(P_gg(:,:,m(j+1)==unique(m)))*V(:,j+1)-c(t(j))*dt-rho(k)*dt],[],2);
end

% Store initial value to measure convergence and perform second
% initialization of secant method:
V_rho(k) = V(g_i,1); rho = [rho 0.9]; k = 2;

V = NaN(length(g),length(t)); V_I = NaN(length(g),length(t));
[V(:,end),V_I(:,end)] = max([g'*R-t_i*rho(k) (1-g)'*R-t_i*rho(k)],[],2);
for j = (length(t)-1):-1:1
    [V(:,j),V_I(:,j)] = max([g'*R-t_i*rho(k) (1-g)'*R-t_i*rho(k) squeeze(P_gg(:,:,m(j+1)==unique(m)))*V(:,j+1)-c(t(j))*dt-rho(k)*dt],[],2);
end
V_rho(k) = V(g_i,1);


% Continue interating using secant method until initial value has
% sufficiently converged:
while abs(V_rho(k)) > tol

    % Update new reward rate using secant method:
    k = k+1; rho(k) = rho(k-1)-V_rho(k-1)*(rho(k-1)-rho(k-2))/(V_rho(k-1)-V_rho(k-2));

    V = NaN(length(g),length(t)); V_I = NaN(length(g),length(t));
    [V(:,end),V_I(:,end)] = max([g'*R-t_i*rho(k) (1-g)'*R-t_i*rho(k)],[],2);
    for j = (length(t)-1):-1:1
        [V(:,j),V_I(:,j)] = max([g'*R-t_i*rho(k) (1-g)'*R-t_i*rho(k) squeeze(P_gg(:,:,m(j+1)==unique(m)))*V(:,j+1)-c(t(j))*dt-rho(k)*dt],[],2);
    end
    V_rho(k) = V(g_i,1);
end
% Return reward rate for converged thresholds:
rho = rho(end);

% Pre-allocate normative thresholds in likelihood space:
g_theta = NaN(1,length(t));

% Construct normative thresholds in likelihood space based off maximal
% index V_I:
for i = length(t):-1:1
    if ~isempty(find(V_I(:,i)==1,1))
        if sum(V_I(:,i)==1)==length(g)
            g_theta(i) = 0.5;
        else
            g_theta(i) = g(find(V_I(:,i)==1,1));
        end
    else
        g_theta(i) = 1;
    end
    if g_theta(i) == 0.5
        g_theta((i+1):end) = 0.5;
    end
end

% Convert normative thresholds to LLR space:
theta = log(g_theta./(1-g_theta));
back to top