Make it possible to use different plotting backends. In particular, add support for Bokeh in addition to Matplotlib.

Hi! I would like to use bayespy in my diploma thesis. It seems more advanced than other libraries available. On the other hand, I spent most of this day reading the documentation trying to figure out how to make a basic network structure according to a given DAG. The closest thing I found was the _add_child() method on Nodes, which is not a public method. Is it okay to use it, to add arbitrary children? As I understand it now, I can't add arbitrary number of parents/children to a node in it's constructor, the number of parents it limited by the number of distribution parameters. I also found the pomegranate library which also implements Bayesian networks and has nice examples for this: https://github.com/jmschrei/pomegranate/blob/master/tutorials/Tutorial_4_Bayesian_Networks.ipynb https://github.com/jmschrei/pomegranate/blob/master/examples/bayesnet_asia.py but pomegranate currently does not implement gaussian nodes, which I'd like to use. Would it be possible to add similar demos to bayespy? Or am I just missing some theory?

Also, I think I'm not the only person who'd appreciate this: http://stackoverflow.com/questions/28431350/create-bayesian-network-and-learn-parameters-with-python3-x

This causes demos to fail!

In [8]: from bayespy import  nodes

In [9]: nodes.Gaussian
Out[9]: bayespy.inference.vmp.nodes.gaussian.Gaussian

In [10]: from bayespy.nodes import  Gaussian
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-10-46b8420f9473> in <module>()
----> 1 from bayespy.nodes import  Gaussian

ImportError: No module named 'bayespy.nodes'

In [11]: import sys

In [12]: sys.ver
sys.version       sys.version_info  

In [12]: sys.version
Out[12]: '3.3.4 |Anaconda 1.9.2 (64-bit)| (default, Feb 10 2014, 17:53:28) \n[GCC 4.1.2 20080704 (Red Hat 4.1.2-54)]'

I guess this should be something trivial?

Running the hmm example in the demos folder, I get the following plotting error at the end of the hmm section:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/usr/local/lib/python3.4/site-packages/matplotlib/colors.py in to_rgba(self, arg, alpha)
    367                     raise ValueError(
--> 368                             'length of rgba sequence should be either 3 or 4')
    369             else:

ValueError: length of rgba sequence should be either 3 or 4

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
/usr/local/lib/python3.4/site-packages/matplotlib/colors.py in to_rgba_array(self, c, alpha)
    398             # Single value? Put it in an array with a single row.
--> 399             return np.array([self.to_rgba(c, alpha)], dtype=np.float)
    400         except ValueError:

/usr/local/lib/python3.4/site-packages/matplotlib/colors.py in to_rgba(self, arg, alpha)
    375             raise ValueError(
--> 376                 'to_rgba: Invalid rgba arg "%s"\n%s' % (str(arg), exc))
    377 

ValueError: to_rgba: Invalid rgba arg "[[  0.00000000e+00   1.00000000e+00   0.00000000e+00]
 [  9.31853843e-04   9.97512196e-01   1.55594968e-03]
...
 [  2.58438900e-03   9.83060709e-03   9.87585004e-01]]"
length of rgba sequence should be either 3 or 4

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
/usr/local/lib/python3.4/site-packages/matplotlib/colors.py in to_rgba(self, arg, alpha)
    364                         raise ValueError(
--> 365                             'number in rbg sequence outside 0-1 range')
    366                 else:

ValueError: number in rbg sequence outside 0-1 range

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-4-157c9bda2cd6> in <module>()
----> 1 run()

<ipython-input-3-c363bc58065a> in run(N, maxiter, seed, std, plot)
     51     colors = Q_hmm['Y'].parents[0]._message_to_child()[0]
     52     plt.plot(y[:,0], y[:,1], 'k-', zorder=-10)
---> 53     plt.scatter(y[:,0], y[:,1], c=colors, s=40)

/usr/local/lib/python3.4/site-packages/matplotlib/pyplot.py in scatter(x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, hold, **kwargs)
   3198         ret = ax.scatter(x, y, s=s, c=c, marker=marker, cmap=cmap, norm=norm,
   3199                          vmin=vmin, vmax=vmax, alpha=alpha,
-> 3200                          linewidths=linewidths, verts=verts, **kwargs)
   3201         draw_if_interactive()
   3202     finally:

/usr/local/lib/python3.4/site-packages/matplotlib/axes/_axes.py in scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, **kwargs)
   3622                 colors = None  # use cmap, norm after collection is created
   3623             else:
-> 3624                 colors = mcolors.colorConverter.to_rgba_array(c, alpha)
   3625 
   3626         faceted = kwargs.pop('faceted', None)

/usr/local/lib/python3.4/site-packages/matplotlib/colors.py in to_rgba_array(self, c, alpha)
    420             result = np.zeros((nc, 4), dtype=np.float)
    421             for i, cc in enumerate(c):
--> 422                 result[i] = self.to_rgba(cc, alpha)
    423             return result
    424 

/usr/local/lib/python3.4/site-packages/matplotlib/colors.py in to_rgba(self, arg, alpha)
    374         except (TypeError, ValueError) as exc:
    375             raise ValueError(
--> 376                 'to_rgba: Invalid rgba arg "%s"\n%s' % (str(arg), exc))
    377 
    378     def to_rgba_array(self, c, alpha=None):

ValueError: to_rgba: Invalid rgba arg "[ 0.  1.  0.]"
number in rbg sequence outside 0-1 range

Running on OS X 10.10, Python 3.4 from Homebrew.

Hi,
Thanks for the great package. I am trying to modify your LDA model code, by using gradient methods (model.optimize()) where I collapse out p_topic and p_word during inference but I am hitting an error.

Basically, in the lda.py example you provided, I changed:
Line 126: Q.update(repeat=300) to the below
126 Q.update(repeat=1)
127 Q.optimize('topics', maxiter=300, collapsed=['p_topic', 'p_word']) # Collapsing p_topic and p_word as we only need individual topic assignments to each word

However I get the following error: "NotImplementedError: Standard gradient not yet implemented for CategoricalDistribution". Is it possible to use BayesPy's advanced methods that support collapsing out variables for models like LDA?

See call stack:

NotImplementedError                       Traceback (most recent call last)
<ipython-input-2-13b3b051706b> in <module>()
----> 1 Q = run(n_documents=100, n_topics=2, n_vocabulary=2, seed=1)

/home/vasant/bayespyexp/lda.py in run(n_documents, n_topics, n_vocabulary, n_words, stochastic, maxiter, seed)
    125 
    126         Q.update(repeat=1)
--> 127         Q.optimize('topics', maxiter=300, collapsed=['p_topic', 'p_word'])
    128 
    129     else:

/home/vasant/.virtualenvs/mypython3/lib/python3.4/site-packages/bayespy/inference/vmp/vmp.py in optimize(self, maxiter, verbose, method, riemannian, collapsed, tol, *nodes)
    451             rg = self.get_gradients(*nodes, euclidian=False)
    452         else:
--> 453             (rg, g) = self.get_gradients(*nodes, euclidian=True)
    454 
    455         if riemannian:

/home/vasant/.virtualenvs/mypython3/lib/python3.4/site-packages/bayespy/inference/vmp/vmp.py in get_gradients(self, euclidian, *nodes)
    370         if euclidian:
    371             g = [self[node].get_gradient(rg_x)
--> 372                  for (node, rg_x) in zip(nodes, rg)]
    373             return (rg, g)
    374         else:

/home/vasant/.virtualenvs/mypython3/lib/python3.4/site-packages/bayespy/inference/vmp/vmp.py in <listcomp>(.0)
    370         if euclidian:
    371             g = [self[node].get_gradient(rg_x)
--> 372                  for (node, rg_x) in zip(nodes, rg)]
    373             return (rg, g)
    374         else:

/home/vasant/.virtualenvs/mypython3/lib/python3.4/site-packages/bayespy/inference/vmp/nodes/expfamily.py in get_gradient(self, rg)
    276         can provide it to this function to avoid re-computing it."""
    277 
--> 278         g = self._distribution.compute_gradient(rg, self.u, self.phi)
    279         for i in range(len(g)):
    280             g[i] /= self.annealing

/home/vasant/.virtualenvs/mypython3/lib/python3.4/site-packages/bayespy/inference/vmp/nodes/expfamily.py in compute_gradient(self, g, u, phi)
     65 
     66         raise NotImplementedError("Standard gradient not yet implemented for %s"
---> 67                                   % (self.__class__.__name__))
     68 
     69 

NotImplementedError: Standard gradient not yet implemented for CategoricalDistribution

Noticed in #34 (comment).

Plate axes are mapped incorrectly in Mixture node. To reproduce, run:

from bayespy.nodes import Categorical, Dirichlet, Mixture
t1 = Categorical([0.2, 0.8])
t2 = Categorical([1/3, 1/3, 1/3])
b = Dirichlet(np.ones((2, 3, 4)))
y = Mixture(t1, Mixture, t2, Categorical, b)
y.observe(3)
b.update()

This will complain about shapes.

Another way to see the problem:

from bayespy.nodes import Categorical, Dirichlet, Mixture
from bayespy.inference import VB

t = Categorical([0.3, 0.7])
b = Dirichlet([[[20, 5], [20, 5]], [[20,5], [5,20]]])
y = Mixture(1, Mixture, t, Categorical, b)
Q = VB(y, b, t)
y.observe(1)
b.update()
print("Incorrect:")
print(b._message_from_children()[0])
b.show()

t = Categorical([0.3, 0.7])
b = Dirichlet([[20,5], [5,20]])
y = Mixture(t, Categorical, b)
Q = VB(y, b, t)
y.observe(1)
b.update()
print("Correct:")
print(b._message_from_children()[0])
b.show()

In the "incorrect" answer, the last 2x2 sub-arrays should be identical to those of the "correct" answer. However, they are not equal.

Currently BayesPy requires Python 3.2, or at least it is tested only on that. Check that it works on 3.3 and add relevant tests.

I wanted to modify the example HMM to additionally estimate the mean and variance of the observations. Oddly, I find that there's rapid convergence to an uninformative model. Any thoughts?

Here's my code (which follows the data generation from the example):

K2 = 10 # >> K
D=2

from bayespy.nodes import Dirichlet
a0 = Dirichlet(1e-3*np.ones(K2))
A = Dirichlet(1e-3*np.ones((K2,K2)))
from bayespy.nodes import CategoricalMarkovChain
Z = CategoricalMarkovChain(a0, A, states=N)

from bayespy.nodes import Gaussian, Wishart
mu_est = Gaussian(np.zeros(D), 1e-5*np.identity(D),plates=(K2,))
Lambda_est = Wishart(D, 1e-5*np.identity(D),plates=())

from bayespy.nodes import Mixture
Y = Mixture(Z, Gaussian, mu_est, Lambda_est)
Y.observe(y)

from bayespy.inference import VB
Q = VB(Y, mu_est, Lambda_est, Z, A, a0)
Q.update(repeat=1000)

Hi, as I try to import bayespy I got this error message:

$ python2.7
Python 2.7.9 (default, Apr  2 2015, 15:33:21) 
[GCC 4.9.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import bayespy
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python2.7/dist-packages/bayespy/__init__.py", line 8, in <module>
    from . import utils
  File "/usr/local/lib/python2.7/dist-packages/bayespy/utils/__init__.py", line 8, in <module>
    from . import misc
  File "/usr/local/lib/python2.7/dist-packages/bayespy/utils/misc.py", line 49
    def parse_command_line_arguments(mandatory_args, *optional_args_list, argv=None):
                                                                             ^
SyntaxError: invalid syntax
>>> 

Give explanation about nested mixtures. It may be unintuitive that the outer mixtures may mix over plate axes that existed in inner mixtures. There is a simple wrapper MultiMixture which avoids this.

Installation from PyPI fails. The command

pip install bayespy

returns an error "no module named numpy" when installing h5py. This happens because either h5py fails to provide proper build requirements or pip fails to solve the order in which dependencies should be installed. See:
h5py/h5py#535
pypa/pip#2478
Suggested solution: drop h5py requirement and make those features extra features. But this is a bit hackish workaround.

Manually installing numpy or h5py before bayespy would solve the issue, but requiring users to manually install some dependencies is not satisfactory because pip should handle that.

The current Beta->Dirichlet converter is confusing because in Beta(a, b), b is the prior sample count for category 1 in Dirichlet but prior sample count for FALSE in Beta. Thus, consider mapping so that Beta([a, b]) is equivalent to Dirichlet([b, a]). This way, Categorical(p) and Bernoulli(p) would be equivalent regardless p is Beta or Dirichlet. See #34.

In any case, add a warning message about the conversion so that users pay extra attention in that case. And so that the warning can be easily silenced with simplefilter.

Currently, the support for discrete directed acyclic graphs (DAGs) is rather poor. There are two options:

• One can construct some simple set of DAGs with Categorical and Mixture nodes, but the posterior approximation is mean-field anyway.
• One can "code" a set of N discrete variables with D possible states into a single discrete variable with D^N states and use Categorical for this node. The inference would be "exact" as the set of discrete variables are not factorized. However, the coding makes things extremely difficult in practice and the inference would not utilize any efficient algorithms (e.g., Junction tree algorithm), thus it would not scale to larger DAGs.

The idea is to implement a node which corresponds to a discrete DAG. It is created by providing conditional probability tables. The node would utilize the Junction tree algorithm for inference. Also, it should be easy to get a single variable or a subset of variables from the node.

X = CategoricalGraph(
    {
        ('rain',): [0.8, 0.2],
        ('sprinkler', 'rain'): [[0.6, 0.4], [0.99, 0.01]],
        ('grass wet', 'sprinkler', 'rain'): [ [[1.0, 0.0], [0.2, 0.8]],
                                              [[0.1, 0.9], [0.01, 0.99]] ],
    }
)

or a bit more complex syntax which allows plates:

X = CategoricalGraph(
    {
        'rain': [0.8, 0.2],
        'sprinkler': dict(table=[[0.6, 0.4], [0.99, 0.01]],
                          given=['rain'],
                          plates=['lawns']),
        'grass wet': dict(table=[ [[1.0, 0.0], [0.2, 0.8]],
                                  [[0.1, 0.9], [0.01, 0.99]] ],
                          given=['sprinkler', 'rain'],
                          plates=['lawns'])
    },
    plates={
        'lawns': 10,
    },
    # If needed, one could explicitly give a mapping from the graph plates
    # to plate axes of other BayesPy nodes:
    plates_axes={
        'lawns': 0,
    }, 
)

Perhaps, the plates could also be specified similarly to other BayesPy nodes, that is, tuples of integers. But then it might be convenient to require that each plate axis is unique meaning that the same plate axis must have the same value throughout the graph. The conditional probability tables can be Dirichlet nodes. One could get a single variable from the DAG as:

rain = X.get_variable('rain')

Or a subset of variables as:

rain_sprinkler = X.get_variables('rain', 'sprinkler')

Observations are given directly to CategoricalGraph:

X.observe('grass wet', 1)

One could maybe even make "probabilistic" observations:

X.observe_probabilities('grass wet', [0.1, 0.9])

which would fix the marginal posterior probabilities.

(Update this sketching as the idea develops.)

First of all, I love this package, the documentation, and especially your (@jluttine) enthusiasm and willingness to help people learn how to use it. You are a huge credit to the community, so thanks!

Second, I have a question. I'm trying to implement the "Mixed Number Subtraction" Bayesian network example discussed in this this paper and in Chapter 11 of the book "Bayesian Networks for Educational Assessment." The basic structure is that there's a 'proficiency model' of latent variables THETA, and an 'evidence model' of measurable/observable variables X. The full problem looks like this:

The latent variables (Skills) are modeled as Bernoulli with Beta priors. Each observable variable (Items) has one or more latent variables that maps to it such for a given observable variable X, if all of the THETA values that map to it have a value of 1, then X has a prior of PI=Beta(20,5). Otherwise, X has a prior PI=Beta(5,20).

The full plate notation looks like this:

I'm omitting the Evidence Model (s) from my implementation until I can get just the i and j plates to work. The delta is a deterministic indicator variable that says, for an item X, whether subject i has all of the prerequisite skills - the prior on X changes if the answer to this (encoded by delta) is true or false. I'm also just setting i=1 to have a single subject.

For a single variable THETA an 10 items I can get the following to work:

lambda1 = Beta([[20,5]])
# using a Categorical here instead of Bernoulli because 
# otherwise it won't work in the Gate below - I get a NoConverterError
theta1 = Categorical(lambda1) 
pi = Beta([[20,5], [5,20]], plates=(10,2)) #two possible beta priors for each for X_1 through X_10
delta = Gate(theta1,pi, gated_plate=-1) #which one is used is based on the value of theta1
X = Binomial(1,delta)
X.observe([1,1,1,1,0,0,0,0,0,0]) #10 separate observations, a single one for each X_i

This seems to work great - i.e. the condition when X depends on the value of only a single THETA. However, I also have the situation in which the prior for theta2 depends on the value of theta1, and the prior for X depends on the value of BOTH variables. What I have so far is:

lambda1 = Dirichlet([20,5]) 
lambda2 = Dirichlet([[5,20],[20,5]])
theta1 = Categorical(lambda1)
theta2 = Categorical(Gate(theta1,lambda2)) #the prior for theta2 depends on the value of theta1

pi1 = Beta([[[5,20], [20,5]]],plates=(10,2)) #two possible beta priors for each question
pi2 = Beta([[[5,20], [20,5]]],plates=(10,2))
delta1 = Gate(theta1,pi1, gated_plate=-1)
delta2 = Gate(theta2,pi2, gated_plate=-1)

X1 = Binomial(1,delta1, plates=(10,))
X2 = Binomial(1,delta2, plates=(10,))

X1.observe([1,1,1,1,1,1,1,1,1,1])
X2.observe([0,0,0,0,0,0,0,0,0,0])

I suspect there's a more compact way to write all the pi1, pi2, etc, but this seems to work. Now what I need to add is a pi3/delta3/X3 where the prior depends on the value of BOTH theta1 and theta2. I suspect the result has to do with a clever use of gates and mixtures like in issue #23, but I don't really understand how you use mixtures with non-Gaussian nodes.

I would really appreciate any guidance! I have this full model implemented in PyMC2 but it's super slow and uses SO much memory. I would really love to switch to variational Bayes using this package.

Edit: Once I get this working I'd be very happy to contribute it as an example in the repo if you're interested.

Edit 2: Fixed an error in the code

initialize_from_value should check that the shape of the array used in initialization is correct. To reproduce the bug, run:

from bayespy.nodes import Gamma
alpha = Gamma(1, 1, plates=(3, 1))
alpha.initialize_from_value(np.ones((1, 3)))
alpha._message_to_child()

This bug was noticed in #36.

Thanks for continued help!

I'm trying to construct an HMM with Poisson observations. Following the approach which works for the Gaussian case, I use a Gamma prior. At the end, I get the error:

Exception: Dimensionality of the observations incorrect.
Shape of input: (300, 3)
Expected shape: (300, 10, 3)
Check plates.

Here's my test case.

import numpy as np
y1 = np.random.poisson([5,0.1,2],size=(100,3))
y2 = np.random.poisson([0.2,6,5],size=(100,3))
y3 = np.random.poisson([2,3,10],size=(100,3))
y = np.vstack([y1,y2,y3])

from bayespy.nodes import Categorical, Dirichlet, Mixture, Poisson, Gamma

beta = 100
alpha = 100
K = 10 # Number of states
(NSamples,D) = y.shape

rate_parameter = Gamma(alpha*np.ones(D), beta*np.ones(D), plates=(K,D))
alpha = Dirichlet(1e-3*np.ones(K), name='alpha')

Z = Categorical(alpha, plates=(NSamples, K, D))
Y = Mixture(Z, Poisson, rate_parameter,cluster_plate=-2)
Y.observe(y)

Thanks again for your rapid replies. I've been trying to trace this issue through the code. In the case where the observation associated with a latent node is Gaussian, we can form an HMM as follows:

# Gaussian HMM
initial_state = Dirichlet(1e-3*np.ones(K))
transmat = Dirichlet(1e-3*np.ones((K,K)))

mu_est = Gaussian(np.zeros(D), 1e-5*np.identity(D),plates=(K,))
Lambda_est = Wishart(D, 1e-5*np.identity(D),plates=(K,))
Y = []
Z = []
for i in range(len(TrainingData[0])):
    [sequence_length, ndims]=TrainingData[0][i].shape
    Z.append(CategoricalMarkovChain(initial_state, transmat, states=sequence_length))
    Y.append(Mixture(Z[i], Gaussian, mu_est, Lambda_est))
    Y[i].observe(TrainingData[0][i])

nodes = Y + [rate_parameter] + Z + [transmat, initial_state]

Q = VB(*nodes)

This works because the Gaussian is implicitly a vector observation distribution. In the case of Poisson if I want to add independent observations (with different densities) to the node, it seems that I should be able to do it like this:

rate_prior = Gamma(alpha, beta, plates=((D,K)))

Y = []
Z = []
for i in range(len(TrainingData[0])):
    [sequence_length, D]=TrainingData[0][i].shape
    Z.append(CategoricalMarkovChain(initial_state, transmat, states=sequence_length))
    Y.append(Mixture(Z[i], Poisson, rate_prior))
    Y[i].observe(TrainingData[0][i])

nodes2 = Y + [rate_parameter] + Z + [transmat, initial_state]
Q2 = VB(*nodes2)

However, this gives me a plate mismatch error: "ValueError: Shapes ([100], (40,)) do not broadcast." (D=100, sequence_length=40). The code works if the "CategoricalMarkovChain" is replaced with a "Categorical".

    Z.append(Categorical(initial_state, plates=(sequence_length, 1)))

I have a feeling that this is easily fixable, but I've been struggling to understand why the plates in the Categorical case are (sequence_length,1), but in the CategoricalMarkovChain case are (sequence_length,).

thanks!

If there is no converter available for a node, print it's name and it's index in the node argument list? See #34 (comment).

Hi,

Congrats on your project. I really like how you organized it.
I am a relative beginner to belief and bayesian networks but I found the documentation and the examples to be relatively easy to understand.
However, I got stuck trying to create a node that might even be a simple one: create a mixture distribution of parents.

I have a set of nodes that result from linear combination of gaussian nodes. Then, I created a new node that would be a simple mixture of the that set of nodes. At first, the Mixture node seemed an obvious choice but the examples in the documentation assume generic Gaussians (or other distributions) and could not understand how to use it when you already have parent nodes.

I also saw from #23 that you can combine several parent nodes in a Mixture node, but could not figure out how to combine to achieve my goal. I could have a the initial argument of the Mixture be the distribution of the clusters, but could not understand how to set the parents.

I sketched and attached a plate notation graphic model figure to try to help explain my goal.

In the figure I have an extra observed node that is a simple sum of the parents. I am trying to figure out how to implement the T^s nodes. I am still exploring different models, and understand how fast it might converge with my data.

Could you provide an example how to do this, if it is possible? Or point me to a point in the documentation that explains this example?
Let me know if I need to make my question more clear and thank you for your time and help!

Add support for non-conjugate nodes and functions. There are a few possibilities of which some or all could be implemented:

• specific variational approximation nodes for some functions as for instance in variational logistic regression
• automatic differentiation variational inference for arbitrary continuous variables and mappings
• black box variational inference for practically anything

Use this issue to sketch the implementation ideas.

When I import bayespy, I got a syntax error:

import bayespy
File "/usr/local/lib/python2.7/dist-packages/bayespy/utils/misc.py", line 452
def sum_multiply(*args, axis=None, sumaxis=True, keepdims=False):
^
SyntaxError: invalid syntax

OS: Ubuntu 12.04 LTS
Python: 2.7
bayespy version: 0.2.3

Imagine that I have observed several different length sequences from a HMM.

At the moment I'm doing the below, which seems quite somewhat annoying in the case in which thousands of sequences are observed. Am I missing something?

Z = CategoricalMarkovChain(a0, A, states=N)
Y = Mixture(Z, Gaussian, mu_est, Lambda_est)
Y.observe(y)

Z1 = CategoricalMarkovChain(a0, A, states=N1)
Y1 = Mixture(Z1, Gaussian, mu_est, Lambda_est)
Y1.observe(y1)

Q = VB(Y, Y1, Z, Z1, mu_est, Lambda_est, A, a0)

Thank you for your work on bayespy, I have found it an extremely useful resource. I am trying to speed up inference for a hierarchical state space model with a constant A (specifically an identity matrix).
I saw that this case was not yet implemented (and noted with a "TODO" in the docs and code) for the parameter expansion by rotation, but was wondering if there was a quick way for me to implement this simplest of cases:

A = np.tile(np.identity(D),(I,1,1,1)) 

where I is simply the number of independent units I observe (and thus y is M x I x N). D is the number of latent dimensions as in your example.
Without the transformation/rotations, I can set everything up (which is thanks to your documentation), but the inference engine proceeds very slowly (as you write).
Thank you for any help you could provide.

Find a maximum likelihood for the concentration parameter in the Dirichlet distribution.

See:

• Thomas P. Minka: "Estimating a Dirichlet distribution"
  http://research.microsoft.com/en-us/um/people/minka/papers/dirichlet/minka-dirichlet.pdf
• Jonathan Huang: "Maximum Likelihood Estimation of Dirichlet Distribution Parameters"
  http://jonathan-huang.org/research/dirichlet/dirichlet.pdf

Write examples as Jupyter notebooks. This makes it easier for readers to try the example themselves.

To use notebooks as a part of a Sphinx-generated documentation, use nbsphinx extension for Sphinx: http://nbsphinx.readthedocs.org/en/latest/

Also, add a direct link to mybinder.org in each notebook so interactive playing with the notebook is even easier. Preferably, also a link to download the notebook itself if Sphinx doesn't generate a nice one automatically.

In order to use these notebooks as a part of doctests, use nosebook: https://github.com/bollwyvl/nosebook. That will ensure that the examples are working.

Not sure why, but it seems one of the code examples raises an error. From the docstring in bernoulli.py...

from bayespy.nodes import Bernoulli, Beta
p = Beta([1e-3, 1e-3])
z = Bernoulli(p, plates=(10,))
z.observe([0, 1, 1, 1, 0, 1, 1, 1, 0, 1])
p.update()
import bayespy.plot as bpplt
import numpy as np
bpplt.pdf(p, np.linspace(0, 1, num=100))

Running this raises the following...

Traceback (most recent call last):
  File "bayespy_example.py", line 8, in <module>
    bpplt.pdf(p, np.linspace(0, 1, num=100))
  File "/path/to/module/bayespy/bayespy/plot.py", line 190, in pdf
    lpdf = Z.logpdf(x)
  File "/path/to/module/bayespy/bayespy/inference/vmp/nodes/expfamily.py", line 500, in logpdf
    (u, f) = self._distribution.compute_fixed_moments_and_f(X, mask=mask)
  File "/path/to/module/bayespy/bayespy/inference/vmp/nodes/beta.py", line 85, in compute_fixed_moments_and_f
    return super().compute_fixed_moments_and_f(p, mask=mask)
  File "/path/to/module/bayespy/bayespy/inference/vmp/nodes/dirichlet.py", line 152, in compute_fixed_moments_and_f
    p = np.asanyarray(p)
UnboundLocalError: local variable 'p' referenced before assignment

Bernoulli mixture model tutorial fails on Python 2.7.10 due to syntax error in "bayespy/utils/misc.py".

Namely, line 49:

def parse_command_line_arguments(mandatory_args, *optional_args_list, argv=None):

Perhaps you meant:

def parse_command_line_arguments(mandatory_args, argv=None, *optional_args_list):

?

This looks like great work. Any thoughts on potentially joining up with PYMC?

Perhaps that could be synergistic-ally beneficial for sets of code, while providing a unified framework for inference in python.

Here is a thread describing a similar idea: https://groups.google.com/forum/#!topic/pymc/Acy9ZdiBRMY

Some code here: https://github.com/kadeng/pypgmc/

I would like to contribute Discrete HMM example (borrowing the wiki example of rainy-sunny).

I was thinking of priors as follows:

theta_start_state ~ beta(20,10)
theta_transition_rainy ~beta(8,2)
theta_transition_sunny ~beta(2,8)
theta_emission_rainy ~ Dirichlet(3,4,3)
theta_emission_sunny ~ Dirichlet(10,6,4)

The following is the code I came up with

# Two states- rainy and sunny
K = 2

# Number of emissions -  walk, shop, clean
V = 3

# Number of observations
N = 100

# Prior on start state
alpha = Dirichlet(np.array([20,10]), name='alpha')

1. I am not clear how to specify A for this trivial example. Should it be defined in terms of plates (as the hmm example you have provided does) or should it be just entered as numpy arrays as I have done above for alpha? Similar question for emission matrix B.
2. How should Y be defined as per the model I stated above? I guess it should be a categorical?

Would be really helpful to add this trivial HMM! Since the mailing list doesn't have many subscribers yet, I decided to pose this to the issue queue directly.

Add more details and examples to Mixture node docstring. See #34

As pointed out in #20, pdf returns incorrect value for observed nodes. Check this. At least for discrete variables, it should work and return 0 or 1. For continuous variables, it might be ok to raise an exception, or return 0 or inf.

I've been trying to implement the following simple stochastic block model:

but I can't figure out how to construct the Mixture node. This is what I have so far:

B = bp.nodes.Beta((1, 1), plates=(N, K))
pi = bp.nodes.Dirichlet(alpha)
q = bp.nodes.Categorical(pi, plates=(N, ))

R = bp.nodes.Mixture(q, bp.nodes.Bernoulli, B)

Any pointers to the right direction?

Take a look at the following code:

import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
import bayespy as bp


def mmsbm_variational(A, K):
    # number of nodes
    N = len(A)

    # directed connection probabilities among groups
    eta = bp.nodes.Beta((1, 1), plates=(K, K))

    # group assignment probabilities for each node
    pi = bp.nodes.Dirichlet(0.5 * np.ones(K), plates=(N,))

    # group assignment for node i connecting to node j
    q_ij = bp.nodes.Categorical(pi, plates=(N,))
    # group assignment for node j connecting to node i
    q_ji = bp.nodes.Categorical(pi, plates=(N,))

    # adjacency matrix
    R = bp.nodes.MultiMixture(
        (q_ij[:, None], q_ji),
        bp.nodes.Bernoulli,
        eta,
        plates=(N, N)
    )

    R.observe(A)

    pi.initialize_from_random()
    Q = bp.inference.VB(R, q_ij, q_ji, pi, eta)
    Q.update(repeat=100)

    return pi.random()


K = 2
n = 10  # number of nodes in each of the 2 clusters
e = n * (n - 1) / 2 * 0.70 # number of edges in each cluster

# create the two separate random graphs and combine them into a single graph
g1 = nx.gnm_random_graph(n, e)
g2 = nx.gnm_random_graph(n, e)
g2 = nx.relabel_nodes(g2, lambda x: x + n)
f = nx.union(g1, g2)

# get the adjacency matrix corresponding to the graph 
A = np.array(nx.adjacency_matrix(f).todense())

# get the mixed membership assignments and draw the with appropriate colors 
mixed = mmsbm_variational(A, K)
nx.draw(f, node_color=[(mixed[i][0], 0, mixed[i][1]) for i in range(len(f.nodes()))])
plt.show()

First I use the networkx library to create a sample graph and visualize it so we can see what's happening.

The important function in terms of bayespy is mmsbm_variational. Here we define the mixed membership stochastic block model.

With this current code I get the expected results (that is we properly separate the nodes in 2 clusters) only around 1/4 times. In other words for different runs I only get some of the time the expected result even though we converge quite soon every single time.

Another thing to notice is that the log-likelihood is significantly lower when we get the expected result compared to when we don't, which makes sense.

Right now I do only the following initialization pi.initialize_from_random(). Any idea if that is what is causing the randomness?

Second question is to what extend does the order of nodes given to bp.inference.VB makes a difference?

Does bayespy suport DBN(Dynamic Bayesian Network), like Factorial HMM?

I am considering switching the license from GPLv3 to MIT.

"the most important two predictors of success for a software project are the number of users and the number of contributors. Because of the restrictions and subtle legal issues involved with GPL licenses, many for-profit companies will not touch GPL-licensed code, even if they are happy to contribute their changes back to the community. A BSD license, on the other hand, removes these restrictions: Hunter mentions several specific examples of vital industry partnership in the case of matplotlib. He argues that in general, a good BSD-licensed project will, by virtue of opening itself to the contribution of private companies, greatly grow its two greatest assets: its user-base and its developer-base."

For more detailed arguments, see:
https://web.stanford.edu/~vcs/talks/VictoriaStoddenIPSC2010.pdf
www.astrobetter.com/blog/2014/03/10/the-whys-and-hows-of-licensing-scientific-code/

If anyone wants to comment on this, please do so quickly, because I'm planning to switch the license this week.

I am still a noob in BayesPy and am trying to implement this variant of LDA described here:
(https://www.cs.cmu.edu/~dbamman/pubs/pdf/bamman+oconnor+smith.acl13.pdf) in Figure 2

Just as in vanilla LDA, we gate the words and their topics based on their topic assignments to words, in this model, words and the corresponding topics are gated but based on assignments to a different variable:entities (that is all words for a given entity need to be drawn from the same topic). Since the number of entities is not the same as the number of words, is it possible to use BayesPy to achieve such a custom gating. When I tried to do this gating, I encountered a broadcasting error (which indicates that I am probably not using this right).

For example below in cell 11: h represents a categorical variable that is gated to psi (selects psi_{p}) if entity e is assigned persona p. Each e is assigned a p (as specified by variable entity_assignments)

# coding: utf-8

# In[1]:

import bayespy
import numpy as np
import bayespy.plot as bpplt
from collections import Counter
import itertools


# In[2]:

# Configuration
numTopics = 2 # Number of Topics 
numPersonas = 2 # Number of Persona. A document is a distribution over personas
numRoles = 3 # Number of roles
numWords = 4 # Vocabulary
numDocuments = 10 # Number of documents
numEntities=4 # Number of entities E in each document. Each entity is assigned a single persona drawn from the document specific distribution over personas
lenEntityDesc = 1000


# In[3]:

# Document specific distribution over personas: theta_{d} is document specific
theta = bayespy.nodes.Dirichlet(np.ones(numPersonas), plates=(numDocuments,))


# In[4]:

entity_doc_ids = []
for docid in np.arange(0, numDocuments):
    entity_doc_ids.extend(itertools.repeat(docid, numEntities))


# In[5]:

entity_doc_ids = np.array(entity_doc_ids).reshape(numDocuments, numEntities)


# In[6]:

# Assignments of personas to each entity in a document. This is the p node in the diagram
entity_assignments = bayespy.nodes.Categorical(bayespy.nodes.Gate(entity_doc_ids, theta), plates=(numDocuments, numEntities))


# In[7]:

# psi_{p}{r} is a distribution over topics. 
psi = bayespy.nodes.Dirichlet(np.ones(numTopics), plates=(numPersonas, numRoles))


# In[8]:

psi.random()


# In[9]:

# phi_{k} is topic which is a distribution over words
phi = bayespy.nodes.Dirichlet(np.ones(numWords), plates=(numTopics,))


# In[10]:

# Role assignments to each word. This is the r node in the diagram.
r = bayespy.nodes.Categorical(np.ones(numRoles)/numRoles, plates=(lenEntityDesc*numEntities,numDocuments))


# In[11]:

# Now I need to gate psi based on personas assigned to each entity. Basically if persona assigned to entity is p and we are looking at
# a word w with role r then we need draw from psi_{p}{r}. So first lets gate on personas and then gate later on roles.
h = bayespy.nodes.Gate(entity_assignments, psi, gated_plate=-2)  

 Encountered a broadcasting error here: ValueError: Shapes ((10, 4), (3,)) do not broadcast

ValueError                                Traceback (most recent call last)
/home/vasant/.virtualenvs/mypython3/lib/python3.4/site-packages/bayespy/inference/vmp/nodes/node.py in _total_plates(cls, plates, *parent_plates)
    288             try:
--> 289                 return misc.broadcasted_shape(*parent_plates)
    290             except ValueError:

/home/vasant/.virtualenvs/mypython3/lib/python3.4/site-packages/bayespy/utils/misc.py in broadcasted_shape(*shapes)
    880                 elif a[i] != 1 and a[i] != s:
--> 881                     raise ValueError("Shapes %s do not broadcast" % (shapes,))
    882         S = S + (s,)

ValueError: Shapes ((10, 4), (3,)) do not broadcast

# In[ ]:
Now gate on r to obtain topic assignments
z = bayespy.nodes.Categorical(bayespy.nodes.Gate(r, h))


# In[ ]:

# Words are drawn from phi based on assigned topics
words = bayespy.nodes.Categorical(bayespy.nodes.Gate(z, phi))

words.random().shape

Hello!

First, a disclaimer, bayesian inference is still quite fuzzy to me, so I may have some foundational misunderstandings that render this question moot.
Also, this is my attempt at going beyond the linear regression example in the docs, so I may have stumbled there as well.
That known though, let's carry on! : )

The Big Picture

I believe this is irrelevant to my problem, but it may highlight some potential misunderstandings of mine, so I will include it. Further, this may be a silly thing to attempt with bayespy. If so, please shoo me out the door. : )
I am trying to fit some data, over a small range in multidimensional space, in the most flexible way possible. Because data at infinity is well outside the relevant range, I've chosen a standard polynomial model. For now, I am pre-calculating the value of each term combination to its respective degree, but that is not included below. Below, I have an arbitrary set of input data "after" the term multiplication and exponent...iation [;)] has been executed.
The first axis of the array is the degree of the term, the second axis differentiates the variable subsets, and the third axis is the value of the term after having the subset multiplied together and being raised to the exponent. The complication is mostly to keep the terms straight in my head. : )

The Specifics

So, I only want to infer the beta values in an array of shape (n, c) #(degree of polynomial, number of combinations). From my hazy understanding, I believe this is accomplished by passing the shape tuple to the plates argument. When I call print(b.random()) it does reveal the expected shape.
Further, the numpy einsum version of what I pass to SumMultiply returns the expected results.
But, SumMultiply complains, via dot.py, that the dimension of b is 0, despite b.random() returning the expected shape.
I tried specifying the ndims argument of GaussianARD, which gets it past the error, but that seems to add those dimensions to the specified shape, and the resulting output fails to line up with the dependent data. I will include the code/output of this if asked.

Again, the snippet below is the simplest version I could write to reproduce the error. Thanks for any and all help!

Code

import bayespy as bp
import numpy as np


n = 3 #the degree of the polynomial to be fitted
d = 3 #the number of dimensions in the independent data
s = 30 #the size each dimension of independent data.

c = 2**d #the number of possible subsets, or combinations, given the number of independent dimensions.

x = np.random.normal(0, 1, (n, c, s)) #arbitrary independent data

b = bp.nodes.GaussianARD(0, 1e-3, plates=(n, c))

print(np.einsum('ij,ijk', b.random(), x))

f = bp.nodes.SumMultiply('ij,ijk', b, x)

Output

[ 352.13537754    6.74444895  230.50129002   26.89123883  -29.92484017
  -33.0472705   137.11028844  112.22195238   -1.93933497 -114.03582012
  -47.26830425 -117.52203611  -16.49295875    8.54611721  147.28272235
  227.64725706 -389.12665218  172.95347527  158.56884408  -31.35686066
 -330.38439688  132.03443101 -251.77565183  178.87004007 -201.6372575
  -91.99130081 -243.31541085   43.45470217  260.13981808  103.8843901 ]
Traceback (most recent call last):
  File "spider/py3/tshoot.py", line 17, in <module>
    f = bp.nodes.SumMultiply('ij,ijk', b, x)
  File "/usr/local/lib/python3.5/site-packages/bayespy/inference/vmp/nodes/dot.py", line 210, in __init__
    len(nodes[n].dims[0])))
ValueError: Wrong number of keys (2) for the node number 0 with 0 dimensions

I've 'pip3 install' bayespy in my virtualenv, but saw a warning when I run nosetest:

---------- begin log ---------------
(myvenv) Weis-MacBook-Air:~ weiliu$ pip freeze
bayespy==0.3.5
Cython==0.22.1
h5py==2.5.0
matplotlib==1.4.3
nose==1.3.7
numpy==1.9.2
pyparsing==2.0.3
python-dateutil==2.4.2
pytz==2015.4
scipy==0.15.1
six==1.9.0
(myvenv) Weis-MacBook-Air:~ weiliu$ nosetests bayespy

................................................................................................E..................

ERROR: Failure: RuntimeWarning (numpy.ufunc size changed, may indicate binary incompatibility)

Traceback (most recent call last):
File "/Users/weiliu/myvenv/lib/python3.4/site-packages/nose/failure.py", line 39, in runTest
raise self.exc_val.with_traceback(self.tb)
File "/Users/weiliu/myvenv/lib/python3.4/site-packages/nose/loader.py", line 418, in loadTestsFromName
addr.filename, addr.module)
File "/Users/weiliu/myvenv/lib/python3.4/site-packages/nose/importer.py", line 47, in importFromPath
return self.importFromDir(dir_path, fqname)
File "/Users/weiliu/myvenv/lib/python3.4/site-packages/nose/importer.py", line 94, in importFromDir
mod = load_module(part_fqname, fh, filename, desc)
File "/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/imp.py", line 235, in load_module
return load_source(name, filename, file)
File "/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/imp.py", line 171, in load_source
module = methods.load()
File "", line 1220, in load
File "", line 1200, in _load_unlocked
File "", line 1129, in _exec
File "", line 1471, in exec_module
File "", line 321, in _call_with_frames_removed
File "/Users/weiliu/myvenv/lib/python3.4/site-packages/bayespy/utils/covfunc/tests/test_distance.py", line 16, in
import scipy.spatial.distance as dist
File "/Users/weiliu/myvenv/lib/python3.4/site-packages/scipy/spatial/init.py", line 90, in
from .ckdtree import *
File "init.pxd", line 861, in init scipy.spatial.ckdtree (scipy/spatial/ckdtree.c:21116)
RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility

Ran 115 tests in 14.507s

FAILED (errors=1)
======= end of log ===============

Any insights about this issue? Seems a upstream issue (bumpy binary compatible), though.

Thanks. And, this package is great work!

Hi,

Great work, sparked my interest in probabilistic programming!

I'm mostly interested in Hidden Markov Models and in particular in a variant called Hidden Markov Trees (not to be confused with hidden markov decision trees). Basically, the hidden states no longer form a chain, but a (binary) tree.
I wondered how I could implement this in bayespy.

My first idea was to create several 'CategoricalMarkovChain' objects and link their hidden states together:

# a0 is the initial distribution of chain a
# P is the prior over transition probabilities
A = CategoricalMarkovChain(a0, P, states=N_A) 

# put the distribution over the last state of A as the initial distribution of B
lastA = A[-1] # somehow get the last variable of chain A (**)
B = CategoricalMarkovChain(lastA, P , states=N_B) 

# same for C
C = CategoricalMarkovChain(lastA, P, states=N_C) 

e.g. chains B,C "inherit" their initial state from chain A.
However, I have no clue how to access the the array of hidden states within a single chain (the line marked with (**) ). Is this possible somehow?

Alternatively one could of course build up the markov chains directly from Categorical, Gate, and Dirichlet variables, such that the hidden states exist explicitly as a list and can be linked between chains as required. However, I'm not sure if inference is still efficient this way (lot more messages/variables probably).

Mixture node requires Categorical parent and doesn't accept Bernoulli parent because there is no converter Bernoulli->Categorical. Add this converter.

This may not be appropriate question to ask in here. If you do know somewhere else I can discuss this kind of questions, please tell me. Thanks!

Fix code style according to PEP8.

Although the requirements explicitly state Python 3.2+, just wanted to confirm if you wish to have bayespy work on Python 2.7.x

I tried a pip install on Python 2.7. The setup gave a few warnings like the following:

      File "/home/nipun/miniconda/lib/python2.7/site-packages/bayespy/inference/vmp/nodes/node.py", line 89
        def __init__(self, *parents, dims=None, plates=None, name=""):
                                        ^
    SyntaxError: invalid syntax

The setup completed fine. However, the same set of errors caused issue while running the demo code.

  File "/home/nipun/miniconda/lib/python2.7/site-packages/bayespy/utils/utils.py", line 367
    def sum_multiply(*args, axis=None, sumaxis=True, keepdims=False):
                               ^
SyntaxError: invalid syntax

I guess this is due to the difference in kwargs and args order handling in Python 2.x vs 3.x?

I am trying to implement a simple mixture of Categorical Distributions where
When p = 0 I need [0,0]
and when p = 1 I need [1,1]

This can easily be achieved through a mixture model as follows: Component Distribution when p =0 should be Categorical([[1,0],[1,0]]) and Component Distribution when p = 1 should be Categorical([[0,1],[0,1]])

p is the selector which is Categorical([0.5,0.5])

However this simple example does not work. It generates a sample [0,1] which should never occurr.

I have verified that individual components generate the right samples. Just trying to use Mixture to select each component however does not work.
.

from bayespy.nodes import Categorical, Mixture
from collections import Counter
p = Categorical([0.5,0.5])

X = np.array([[[0.999,0.001],[0.999,0.001]],[[0.001,0.999],[0.001,0.999]]])

z.random()
Output: array([0, 1]) << We should never actually see (0,1) because when we select a component it should only either yield [0,0] or [1,1]

q = Categorical(X[0], plates=(1000,2))  << This works only [0,0] should be returned
array([[0, 0],
       [0, 0],
       [0, 0],
       ..., 
       [0, 0],
       [0, 0],
       [0, 0]])

Similarly X[1] Categorical(X[1], plates=(1000,2)) only returns [1,1]

Support exact inference for (sub-)graphs of categorical variables by using the junction tree algorithm.

BernoulliMoments arrays can be converted to GaussianMoments arrays, and thus BernoulliArray could be used to mask Gaussian arrays by multiplying. Enables, for instance, spike and slab priors.

Pretty new to bayesian graphical models (and your library) so this may well be a n00b question...

I'm looking to model what I guess would be called a markov tree (not hidden, in contrast to one of the other discussions). Each node has a conditional probability table, with the likelihood of being in each state dependent on the state of the parent. In line with the discussion here I can construct a hidden markov tree, but not (as I understand it) then fix the true state of any given node. Right?

If I want to take the most simplistic approach and just use the parent node in the probability array of the child, how can I use a categorical node? Is this what a gate node is for?! (I couldn't find them used anywhere in an example...)

Does that question make any sense...?!

Are there any examples of using BayesPy to do multilevel logistic regression? This seems like it'd be straightforward to do if you could define a variable as an arbitrary function of another, like

p[i] <- 1 / (1 + exp(-z[i]))

but I don't know whether that's possible.