Commit dc6d19ff authored by Tiago Peixoto's avatar Tiago Peixoto

blockmodel: Fix parameter names of weight distributions and add support for Jeffreys prior

parent 477df9b8
Pipeline #265 passed with stage
in 390 minutes and 12 seconds
......@@ -258,6 +258,8 @@ double positive_w_log_P(DT N, double x, double alpha, double beta)
{
if (N == 0)
return 0.;
if (alpha == 0 && beta == 0) // Jeffreys
return lgamma(N) - N * log(x);
return lgamma(N + alpha) - lgamma(alpha) + alpha * log(beta) -
(alpha + N) * log(beta + x);
}
......@@ -269,9 +271,18 @@ double signed_w_log_P(DT N, double x, double v, double m0, double k0, double v0,
{
if (N == 0)
return 0.;
if (v0 == 0 && k0 == 0) // Jeffreys
{
if (N > 1)
return lgamma(N/2.) - log(N)/2 - N * log(v) / 2 - (N / 2.) *
log(M_PI);
else
return -log(pow(x - m0, 2) + v0)/2;
}
auto k_n = k0 + N;
auto nu_n = nu0 + N;
auto v_n = (v0 * nu0 + v + ((N * k0)/(k0 + N)) * pow(m0 - x/N, 2.)) / nu_n;
auto v_n = (v0 * nu0 + v + ((N * k0)/(k0 + N)) * pow(m0 - x/N, 2)) / nu_n;
return lgamma(nu_n / 2.) - lgamma(nu0 / 2.) + (log(k0) - log(k_n)) / 2. +
(nu0 / 2.) * log(nu0 * v0) - (nu_n / 2.) * log(nu_n * v_n) - (N / 2.) *
log(M_PI);
......@@ -283,6 +294,8 @@ double geometric_w_log_P(DT N, double x, double alpha, double beta)
{
if (N == 0)
return 0.;
if (alpha == 0 && beta == .5) // Jeffreys
return lbeta(N, x + .5);
return lbeta(N + alpha, x + beta) - lbeta(alpha, beta);
}
......@@ -297,12 +310,14 @@ double binomial_w_log_P(DT N, double x, size_t n, double alpha, double beta)
// discrete: Poisson
template <class DT>
double poisson_w_log_P(DT N, double x, double r, double theta)
double poisson_w_log_P(DT N, double x, double alpha, double beta)
{
if (N == 0)
return 0.;
return lgamma(x + r) - lgamma(r) - r * log(theta) - (x + r) *
log(N + 1. / theta);
if (alpha == .5 && beta == 0) // Jeffreys
return lgamma(x + .5) - (x + .5) * log(N) - lgamma(.5);
return lgamma(x + alpha) - (x + alpha) * log(N + beta) - lgamma(alpha) +
alpha * log(beta);
}
// ===============
......
......@@ -448,9 +448,10 @@ class BlockState(object):
self.wparams = libcore.Vector_Vector_double()
for i, rt in enumerate(self.rec_types):
ps = Vector_double()
if rt == libinference.rec_type.real_exponential:
defaults = OrderedDict([("r", 1),
("theta", self.recs[i].fa.mean())])
if rt in [libinference.rec_type.real_exponential,
libinference.rec_type.discrete_poisson]:
defaults = OrderedDict([("alpha", 1),
("beta", self.recs[i].fa.mean())])
elif rt == libinference.rec_type.real_normal:
defaults = OrderedDict([("m0", self.recs[i].fa.mean()),
("k0", 1),
......@@ -459,9 +460,6 @@ class BlockState(object):
elif rt == libinference.rec_type.discrete_geometric:
defaults = OrderedDict([("alpha",1),
("beta",1)])
elif rt == libinference.rec_type.discrete_poisson:
defaults = OrderedDict([("r", 1),
("theta", self.recs[i].fa.mean())])
elif rt == libinference.rec_type.discrete_binomial:
defaults = OrderedDict([("N", self.recs[i].fa.max()),
("alpha", 1),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment