An asymmetric pdf with infinite support

When was the last time you needed an asymmetric (skewed) probability density function (pdf) with infinite support? Traditional skewed  distributions like the gamma family suffer from a semi-infinite support, that is, \mathrm{supp}(p) = [0, +\infty). The support, if you are out of the loop, is the set of values x in the domain of p such that p(x) > 0. Why is this inconvenient? well, I will give more details about the specific application later, meanwhile, let’s say that the fact that its derivative is discontinuous at 0 is problematic; even more, I need my function to be at least in C^3, that is, to have at least 3 continuous derivatives!.

This post constitutes a somewhat dirty solution being as unwilling as I am to review any literature in depth (I might be inventing the wheel again, some wheel, but who cares; this is a blog!). Here we go.

Firstly, it is important that p(t) > 0 for all t. I also need that \lim_{t \rightarrow \pm\infty} = 0. One function that satisfies this, at least for t > 0 is the so-called \alpha-function (this name for this particular function seems to be exclusively used in neuroscience). It has the form:

p(t;A, \mu, \tau) = A(t - \mu)\exp(-k (t - \mu)),

defined only for t > 0. k is actually the ratio of the rise/decay rates of the curve; you can think about that as some sort of spread of the distribution. The graph of the \alpha-function looks something like this:

alpha

From now on, the domain of the pdf will be time (but t \in (-\infty, \infty)!). In order to get the other half of the domain filled, I will glue a second \alpha-function but supported in (-\infty, 0]. With the appropriate glue, the resulting function should be supported in the whole real line. How does this “appropriate” “glue” looks like? Well, let’s start by imposing some conditions. Let our two \alpha-functions be: p_1(t; A_1, \mu_1, k_1) and p_2(-t;A_2, \mu_2, k_2) (notice the inverted the time for the second one). The conditions are then:

  1. There should be some \tau such that p_1(\tau) = p_2(\tau)
  2. In order for the connection to be “smooth”, p_1'(\tau) = p_2'(\tau), where ' = \frac{d}{dt}. That is, the derivatives should be the same at that point so that one function follows smoothly after the other.
  3. Looking at the domains of p_1 and p_2, considering that \tau should be always within the domain of both functions, I should have \mu_1 < \tau < \mu_2
  4. Ideally, also the higher order derivatives p_1^{(n)}(\tau) = p_2^{(n)}(\tau) for n at most 3, but I won’t worry about that for the time being.

Now, I will follow a semi-analytical procedure to find \mu_2, \tau and A_2, given k_1, k_2, \mu_1 and A_1 = 1. Effectively, what I am going to do is to keep p_1 still while translating and scaling p_2 until it satisfy the conditions above. Instead of trying to bruteforce my way through the problem, let’s see what kind of information can be extracted from the previous conditions.

From the condition 1, (\tau - \mu_1)\exp(-k_1(\tau - \mu_1)) = A (\mu_2 -\tau)\exp(-k_2( \mu_2-\tau)),we have

A = \frac{\tau - \mu_1}{\mu_2 - \tau}\exp[-(k_1 + k_2)\tau + k_1\mu_1 + k_2\mu_2].

From the condition 2,

A = -\frac{1 - k_1(\tau - \mu_1)}{1 - k_2(\mu_2 - \tau)}\exp[-(k_1 + k_2)\tau + k_1\mu_1 + k_2\mu_2].

Combining those two, I conclude that, for the right values of  \hat{\mu_2} and \hat{\tau}, the rational function

g (\tau,\mu_2) = \frac{\tau - \mu_1}{\mu_2 - \tau} + \frac{\tau - \mu_1}{\mu_2 - \tau},

should be zero (i.e. g(\hat{\tau}, \hat{\mu_2}) = 0). I could keep looking through the next derivatives for more equations to get the exact values, but, as lazy as I am, I switch to numerical mode. Using the fsolve function of Matlab,  let’s solve iteratively the following surrogate function

f(\tau) = g(\tau, x)

for each x \in (\mu_1, M) for some M. As it can be seen in the following picture, some times the function will have no roots (above the zero line)

functions1

that means that that particular value of \mu_2 is not good and we move on. Note the two asymptotes of f at \tau = \mu_2 and  \tau = \mu_2 - \frac{1}{k_2}.

functions2

The value of \mu_2 for which a root is found corresponds to a good glue point. We can now go on and compute A and combine the p_1 and p_2 in a final function p(t) = \Theta(t -\tau)p_1 + \Theta(\tau-t)p_2:

pFInally, it remains to be seen if this function has more derivatives (we might want to add more restrictions). Perhaps, we would also want to normalize it such that \int_{\mathbb{R}} p(t) \mathrm{d}t = 1. Besides, we could try to find the new mean, variance, kurtosis and the like. Those things will have to wait until a future post.

A short word about the application. I wanted a highly smooth synthetic firing rate density for neurons with that prescribed profile (for more information about this particular detail, check the corresponding post… in the future).  This profile is not expected to be seen (sought after?) in data (so no data fitting available), it is just a hypothetical profile that was desired to answer some questions that go beyond the current topic.

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s