# Is there an analytic solution for the density function of this complex random variable?

The process below yields a distribution of "response times" (RT), and I'd like to know if there is an analytic solution to obtain the density function of this distribution.

An RT is recorded at the earliest time ($t_1$) at which a function, $f(t)$, reaches a threshold.

The threshold varies across iterations of the process according to a random uniform distribution with bounds specified by the parameters $P_1$ and $P_2$.

$f(t)$ begins at $t=0$ and has a form whose derivative is defined by the function $g(t)$:

$$g(t) = P_3\left (\frac{P_4}{P_4+P_5}\right )\left (1-\exp\left [-(P_4+P_5)\int_{0}^{t}h(s)ds\right ]\right )$$

where $h(t)$ is defined as:

$$h(t) = P_6\Gamma(t,P_8,n)[1 - \Gamma(t-P_7,P_9,n)]$$

and the gamma function defined as:

$$\Gamma(t,\beta,n) = 1- \left ( e^{-\beta t} \right )\left ( \sum_{n-1}^{j=0} \frac{(\beta t)^{j}}{j!} \right )$$

So, there are in total 10 parameters in the process: $P_1$ through $P_9$, and $n$. I could of course perform a Euler simulation to generate a distribution of RT values for given a particular set of 10 parameter values, but it would be nice to avoid this because I'm hoping to use a Nelder–Mead search algorithm to find the set of 10 parameter values that maximizes the likelihood of a set of observed RT data points, and it really slows down the search if I have to compute estimates of the density function for every candidate set of parameter values.

Pursuant to comments, here is a further description of the process:

The process is an adaptation of the model presented by Smith & Racliff (PDF) which models perception, attention and decision in the context of a speeded detection task. Function $h(t)$ models perception: accumulation of information about a brief visual stimulus. Function $g(t)$ models the transfer of this information to memory via attention, and function $h(t)$ is a linear ballistic decision process in the vein of Brown & Heathcote (PDF). I thought that supplementing an LBA as the decision mechanism in place of Smith & Ratcliff's diffusion process might simplify the derivation of the final density function.

$P_1$ and $P_2$ are the upper and lower bounds of the random uniform threshold. $P_3$ is the asymptote of $g(t)$. $P_4$ is the intensity of the stimulus. $P_5$ is the intensity of the background. $P_6$ is the amplitude of the stimulus. $P_7$ is the duration of the stimulus. $P_8$ is the rate at which the stimulus comes on. $P_9$ is the rate at which the stimulus turns off.

Smith & Ratcliff (2009) note that $n$ might be reasonably fixed to a value of 3, so possibly that can further simplify things. Additionally, while Smith & Ratcliff sought to fit their model to different stimulus intensities, I don't need this particular feature, so I'm not sure if I need parameters $P_4$ and $P_5$ (though I'm not sure how to properly remove them either).

-
For some reason, I have the funny feeling that all the parameters from $P_1$ to $P_9$ aren't completely independent. Mike, you should probably describe what each of those nine parameters represent, it can help with the analysis. –  Guess who it is. Aug 21 '10 at 14:21
@J. Mangaldan: Point taken; I updated the question to include more information about the process and parameters. –  Mike Lawrence Aug 22 '10 at 0:14

I can't find a closed-form, but I can get partway. Note that your $\Gamma$ can be expressed as $1-\Gamma(n,\beta t)/\Gamma(n)$ using the incomplete gamma function; it can also be expressed (in more complicated form) using the WhittakerM function (aka 1F1). Thus so can $h(t)$. What this really means is that $h(t)$ is holonomic and it particular satisfies a 4th order linear ordinary differential equation with polynomial coefficients. The product-integral in $g$'s definition bumps the order by one, and then integration (to get $f$) by another order, so $f(t)$ satisfies an order $6$ LODE. With all those parameters, it is hopeless to solve.

If you have a few relations - like if $P_8$ and $P_9$ are related - then you might have a chance. Basically, unless you can knock down the complexity of what $h$ actually is, it is rather unlikely that the whole process will 'collapse' in the presence of that many free parameters.

Assuming that $n=3$ makes things much simpler. Then $h(t)$ is a simple sum of exponentials. It is however integrable in closed-form (but is a medium-sized mess, so I won't post it here). But that's where things are likely to stop, since exponentials-of-exponentials do not tend to integrate in closed-form very often. Numerically, one ought to be able to generate some very good approximations for evaluating $f(t)$ though. I guess it depends on what you are eventually trying to do with your equations.

-

I can't hope to do better than Jacques's answer, so I'll just offer some little numerical advice on computing your function to be fitted (but you probably already know this): make sure that whatever system you're doing your computations in has support for both the upper incomplete gamma function $\Gamma(a,x)$ and the lower incomplete gamma function $\gamma(a,x)$. With these two, we can represent your $h(t)$ as

$h(t)=\frac1{(n-1)!^2}P_6 \gamma(n,P_8 t)\Gamma(n,P_9(t-P_7))$

The reason for doing things this way is that we mitigate the possibility of any subtractive cancellation happening when your downhill simplex optimization wanders too far out. Which brings me to my second piece of advice: you would do well to find good initial values for your $P_i$ parameters. Nelder-Mead isn't really known for being fast (at least there's no need to assume continuity!), and it also tends to wander around your parameter space unless suitably restricted.

-
Excellent suggestion regarding upper/lower incomplete gammas. –  Jacques Carette Aug 22 '10 at 14:25
I should probably also add: if the system supports regularized versions of both incomplete gamma functions (the versions where the factorial denominators are incorporated), it's better to use those than the vanilla incomplete gammas. –  Guess who it is. Aug 22 '10 at 15:42