[SciPy-User] Joint distributions
William Furnass
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
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
josef...@gmail.com
scipy.stats only has univariate distributions, or to be exact it
calculates it for many points independently.
So the returned array is the pdf for each point separately calculated.
If you want the pdf for the bivariate or multivariate normal
distribution then it's just a few lines,
( I think the bivariate normal is also in matplotlib, in statsmodels ?)
Your fitting problem sounds a bit like what scipy.odr does.
Josef