A function that uses model selection to detect variability in the baseline firing rate of a cell. The algorithm assumes Poisson spiking in the baseline period and compares the posterior on a single rate explaining the data to having two rates.

This is an approximation for the ideal situation of comparing a single rate to the possibility of having more than one rate, i.e. two,three,four, etc. That is, given the vector of baseline rates r, we're saying

DetectVariableBaselineUsingBayesianModelSelection(r, lrange = NULL,
  nsd = 5, nlvals = 1001, bias = 0)

Arguments

r	A vector of length N contain the baseline spike counts for N odors.
lrange	The range of baseline rates to integrate over. If lrange is NULL, will compute this from the data.
nsd	The number of standard deviations around the observed range of rates to integrate over.
nlvals	The number of points along each dimension to use in perform the numerical integration.
bias	The bias term in favour of M1.

Value

A list consisting of

bestModel

1 or 2, indicating the better model.

lpM1

The log likelihood of the single rate model.

lpM2

The log likelihood of the two rate model.

Details

relatve prob. of single rate = p(M_1|r)/(p(M_2|r) + p(M_3|r) + ...) ~ p(M_1|r)/p(M_2|r).

Hence the procedure is theoretically biased towards the single rate model, and the user can supply a counter bias to this in the bias argument.

Comparing the two models requires a prior probability on each:

p(M_1|r) / p(M_2|r) = p(r|M_1)p(M_1)/p(r|M_2)p(M_2)

We actually compare the logarithm of this

log [p(M_1|r)/p(M_2|r)] = log [p(r|M_1)/p(r|M_2)] + log [p(M_1)/p(M_2)]

We assume p(M_1) = p(M_2) so that the default comparison is that of the likelihoods, and the latter term is zero, but the user can provide the logarithm of the prior probability ration in the bias term. The final decision of the model is:

M1 if log [p(r|M_1)/p(r|M_2)] + bias > 0, otherwise M2.

The likelihoods are computed by marginalizing over the rates. For example, for the two rate case, we have

p(r|M_2) = 0.5 * int(l1, lmin, lmax) int(l0, lmin, lmax) p(r|l0,l1,M_2) p(l0)p(l1)

The prior on the rates are assumed uniform, i.e. p(l0) = p(l1) = 1/(lmax - lmin).

The integration is performed numerically by splitting the range (lmin-lmax) into a fixed pitch grid, computing the probability at each grid point, and summing the result scaled by the grid pitch.