function [v,a,T] = EZdiff_4choice(VRT,Pc,MRT)

% Author: Jeffrey M. Spielberg (jspielb2@gmail.com)
% Version: 01.06.15
% 
%   Calculates diffusion model parameters based on Wagenmakers et al.
%   (2007)'s EZ-diffusion model with one change. Specifically, Wagenmakers
%   et al. assume that there is no bias among the response options. In a 
%   two-option task, this would mean that the starting point (z) should be 
%   halfway between the two boundary edges (i.e., z = a/2). We can also use
%   the diffusion model for 4-option tasks if we restructure the output to 
%   be correct vs. incorrect (rather than option A vs. option B). In this 
%   case, assuming that there is no bias among the response options means
%   that there is a 3 to 1 bias against the correct option (i.e., z = a/4).
%   This code calculates the diffusion parameters using this assumption.
% 
%   NOTE, diffusion models in general, and the EZ-diffusion model in
%   particular, have not been extensively tested for use in tasks with 2+
%   response options. Thus, the following code should be used with caution.
%
%   Usage: [v,a,T] = EZdiff_4choice(VRT,Pc,MRT)
% 
%       VRT = variance of reaction times for correct trials
%       Pc  = percent correct
%       MRT = mean reaction time for correct trials
% 
%       v = mean drift rate
%       a = boundary separation
%       T = mean non-decision time
% 
%   Parameters for each condition should be calculated separately. 
% 
% 
% WARNING: This is a beta version. There no known bugs, but only limited
% testing has been perfomed. This software comes with no warranty (even the
% implied warranty of merchantability or fitness for a particular purpose).
% Therefore, USE AT YOUR OWN RISK!!!
%
% Copyleft 2015. Software can be modified and redistributed, but modifed,
% redistributed versions must have the same rights

s      = 0.1;
lambda = (Pc-1)/Pc;
gamma  = (7-(27*lambda))/54;
cgamma = nthroot(((6+(54*(lambda^2))-(28*lambda))/216),2);
alpha  = nthroot((gamma+cgamma),3);
beta   = nthroot((gamma-cgamma),3);
kappa  = alpha+beta-(1/3);

v = nthroot(((-(s^4)*log(kappa))/VRT)*((4*log(kappa)*(kappa^2)-(kappa^4)+1)/(((kappa^2)+1)^2)),4);
a = (-2*(s^2)*log(kappa))/v;
T = MRT-((a/(2*v))*((1-(kappa^2))/(1+(kappa^2))));
