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))));