[SciPy-User] Joint distributions

[William Furnass]

I am wanting to fit a parameterised model to a series of 15
datapoints, with each being a concentration C and time t.  Within the
objective function of the optimisation routine that I'm using for the
model fitting I presently calculate fitness using the Bray Curtis
distance between the data series and the prediction corresponding to a
candidate solution.

I would ideally like to calculate fitness in such a way as to account
for uncertainty in each (C, t).  I think I can achieve this for a
given data series by
 a) modelling each data point using a bivariate Gaussian PDF (with
static variances for both C and t)
 b) calculate a prediction using a small dt
 c) find the highest probability of all points in the prediction
series for each of the 15 bivariate PDFs
 d) sum or average the probabilities to get a measure of the fit of
the real data series to the prediction corresponding to the candidate
solution.

My question is is there an easy way of finding joint probabilities
using scipy.stats?  I thought I could construct a bivariate normal
distribution using

dens = scipy.stats.norm(loc=np.array([t[i], C[i]]),
scale=np.array([t_stdev, C_stdev]))

but

dens.pdf(np.array([5,7]))

returns an array when I thought it should return a scalar probability.

Apologies if the above is not particularly clear or if I'm missing
something obvious here.

Regards,

Will Furnass