{"action":"created","issue":{"url":"https://api.github.com/repos/bayespy/bayespy/issues/52","repository_url":"https://api.github.com/repos/bayespy/bayespy","labels_url":"https://api.github.com/repos/bayespy/bayespy/issues/52/labels{/name}","comments_url":"https://api.github.com/repos/bayespy/bayespy/issues/52/comments","events_url":"https://api.github.com/repos/bayespy/bayespy/issues/52/events","html_url":"https://github.com/bayespy/bayespy/issues/52","id":139235880,"number":52,"title":"How to do proper initialization?","user":{"login":"abojchevski","id":8150330,"avatar_url":"https://avatars.githubusercontent.com/u/8150330?v=3","gravatar_id":"","url":"https://api.github.com/users/abojchevski","html_url":"https://github.com/abojchevski","followers_url":"https://api.github.com/users/abojchevski/followers","following_url":"https://api.github.com/users/abojchevski/following{/other_user}","gists_url":"https://api.github.com/users/abojchevski/gists{/gist_id}","starred_url":"https://api.github.com/users/abojchevski/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/abojchevski/subscriptions","organizations_url":"https://api.github.com/users/abojchevski/orgs","repos_url":"https://api.github.com/users/abojchevski/repos","events_url":"https://api.github.com/users/abojchevski/events{/privacy}","received_events_url":"https://api.github.com/users/abojchevski/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":1,"created_at":"2016-03-08T10:17:22Z","updated_at":"2016-03-10T07:34:41Z","closed_at":null,"body":"Take a look at the following code:\r\n\r\n```\r\nimport networkx as nx\r\nimport numpy as np\r\nimport matplotlib.pyplot as plt\r\nimport bayespy as bp\r\n\r\n\r\ndef mmsbm_variational(A, K):\r\n    # number of nodes\r\n    N = len(A)\r\n\r\n    # directed connection probabilities among groups\r\n    eta = bp.nodes.Beta((1, 1), plates=(K, K))\r\n\r\n    # group assignment probabilities for each node\r\n    pi = bp.nodes.Dirichlet(0.5 * np.ones(K), plates=(N,))\r\n\r\n    # group assignment for node i connecting to node j\r\n    q_ij = bp.nodes.Categorical(pi, plates=(N,))\r\n    # group assignment for node j connecting to node i\r\n    q_ji = bp.nodes.Categorical(pi, plates=(N,))\r\n\r\n    # adjacency matrix\r\n    R = bp.nodes.MultiMixture(\r\n        (q_ij[:, None], q_ji),\r\n        bp.nodes.Bernoulli,\r\n        eta,\r\n        plates=(N, N)\r\n    )\r\n\r\n    R.observe(A)\r\n\r\n    pi.initialize_from_random()\r\n    Q = bp.inference.VB(R, q_ij, q_ji, pi, eta)\r\n    Q.update(repeat=100)\r\n\r\n    return pi.random()\r\n\r\n\r\nK = 2\r\nn = 10  # number of nodes in each of the 2 clusters\r\ne = n * (n - 1) / 2 * 0.70 # number of edges in each cluster\r\n\r\n# create the two separate random graphs and combine them into a single graph\r\ng1 = nx.gnm_random_graph(n, e)\r\ng2 = nx.gnm_random_graph(n, e)\r\ng2 = nx.relabel_nodes(g2, lambda x: x + n)\r\nf = nx.union(g1, g2)\r\n\r\n# get the adjacency matrix corresponding to the graph \r\nA = np.array(nx.adjacency_matrix(f).todense())\r\n\r\n# get the mixed membership assignments and draw the with appropriate colors \r\nmixed = mmsbm_variational(A, K)\r\nnx.draw(f, node_color=[(mixed[i][0], 0, mixed[i][1]) for i in range(len(f.nodes()))])\r\nplt.show()\r\n```\r\n\r\nFirst I use the networkx library to create a sample graph and visualize it so we can see what's happening.\r\n\r\nThe important function in terms of bayespy is `mmsbm_variational`. Here we define the mixed membership stochastic block model.\r\n\r\nWith 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.\r\n\r\nAnother 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.\r\n\r\nRight now I do only the following initialization `pi.initialize_from_random()`. Any idea if that is what is causing the randomness?\r\n\r\nSecond question is to what extend does the order of nodes given to bp.inference.VB makes a difference?\r\n\r\n\r\n\r\n\r\n"},"comment":{"url":"https://api.github.com/repos/bayespy/bayespy/issues/comments/194716542","html_url":"https://github.com/bayespy/bayespy/issues/52#issuecomment-194716542","issue_url":"https://api.github.com/repos/bayespy/bayespy/issues/52","id":194716542,"user":{"login":"jluttine","id":2195834,"avatar_url":"https://avatars.githubusercontent.com/u/2195834?v=3","gravatar_id":"","url":"https://api.github.com/users/jluttine","html_url":"https://github.com/jluttine","followers_url":"https://api.github.com/users/jluttine/followers","following_url":"https://api.github.com/users/jluttine/following{/other_user}","gists_url":"https://api.github.com/users/jluttine/gists{/gist_id}","starred_url":"https://api.github.com/users/jluttine/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/jluttine/subscriptions","organizations_url":"https://api.github.com/users/jluttine/orgs","repos_url":"https://api.github.com/users/jluttine/repos","events_url":"https://api.github.com/users/jluttine/events{/privacy}","received_events_url":"https://api.github.com/users/jluttine/received_events","type":"User","site_admin":false},"created_at":"2016-03-10T07:34:41Z","updated_at":"2016-03-10T07:34:41Z","body":"\r\nHere are some thoughts on your implementation:\r\n\r\n- You are using `plates=(N,)` for `pi`. It might be what you want, but another option is to discard the plates. That way the prior probability to belong to a group depends on how many other nodes are in the same group, and some very small groups may be pruned out. In addition, you could use different `pi` for `q_ij` and `q_ji`. But this was just a thought, it really depends on what you want.\r\n\r\n- About the initialization: reasonable initialization is sometimes a crucial step in VB (and in MCMC too). Some general rules of thumbs I'd say:\r\n\r\n   - use randomness when need to break symmetry\r\n   - use some given values when you want to give a little push to the right direction\r\n   - use these custom initializations typically for higher level hyperparameters or smaller nodes\r\n   - update \"uninitialized\" nodes (well, BayesPy does a simple initialization for all nodes) and/or lower level larger nodes first by giving them first in the node list for `VB` object\r\n   \r\n- Thus, I made the following changes to the initialization:\r\n\r\n   - initialize `eta` with a probability array which aims at grouping mutually connected nodes together\r\n   - use randomness to break symmetry in `q` variables, and update them first\r\n   \r\n- Your MMSBM function return random sample of `pi`. That might be what you wanted if you wanted to use `plates=(N,)` for it. But another option is to return directly the posterior group assignment probabilities of each node by `q_ij.get_moments()[0]`. Also, `q_ji.get_moments()[0]` would give the other group assignment probabilities for each node.\r\n\r\n\r\n```python\r\nimport networkx as nx\r\nimport numpy as np\r\nimport matplotlib.pyplot as plt\r\nimport bayespy as bp\r\n\r\n\r\ndef mmsbm_variational(A, K):\r\n    # number of nodes\r\n    N = len(A)\r\n\r\n    # directed connection probabilities among groups\r\n    eta = bp.nodes.Beta((1, 1), plates=(K, K))\r\n\r\n    # group assignment probabilities for each node\r\n    pi = bp.nodes.Dirichlet(0.5 * np.ones(K))  # perhaps discard the plates?\r\n    #pi = bp.nodes.Dirichlet(0.5 * np.ones(K), plates=(N,))\r\n\r\n    # group assignment for node i connecting to node j\r\n    q_ij = bp.nodes.Categorical(pi, plates=(N,))\r\n    # group assignment for node j connecting to node i\r\n    q_ji = bp.nodes.Categorical(pi, plates=(N,))  # you could use different pi here?\r\n\r\n    # adjacency matrix\r\n    R = bp.nodes.MultiMixture(\r\n        (q_ij[:, None], q_ji),\r\n        bp.nodes.Bernoulli,\r\n        eta,\r\n        plates=(N, N)\r\n    )\r\n\r\n    R.observe(A)\r\n\r\n    eta_init = [[0.9 if i == j else 0.05 for i in range(K)] for j in range(K)]\r\n    eta.initialize_from_value(eta_init)\r\n    q_ij.initialize_from_random()\r\n    q_ji.initialize_from_random()\r\n    Q = bp.inference.VB(R, q_ij, q_ji, eta, pi)\r\n    Q.update(repeat=100)\r\n    \r\n    plt.figure()\r\n    bp.plot.hinton(eta)\r\n\r\n    return q_ij.get_moments()[0]\r\n    # return q_ji.get_moments()[0]  # you might be interested in this too\r\n    #return pi.random()  # were you really interested in `pi`?\r\n\r\n\r\nK = 2\r\nn = 10  # number of nodes in each of the 2 clusters\r\ne = n * (n - 1) / 2 * 0.70 # number of edges in each cluster\r\n\r\n# create the two separate random graphs and combine them into a single graph\r\ng1 = nx.gnm_random_graph(n, e)\r\ng2 = nx.gnm_random_graph(n, e)\r\ng2 = nx.relabel_nodes(g2, lambda x: x + n)\r\nf = nx.union(g1, g2)\r\n\r\n# get the adjacency matrix corresponding to the graph \r\nA = np.array(nx.adjacency_matrix(f).todense())\r\n\r\n# get the mixed membership assignments and draw the with appropriate colors \r\nmixed = mmsbm_variational(A, K)\r\n\r\nplt.figure()\r\nnx.draw(f, node_color=[(mixed[i][0], 0, mixed[i][1]) for i in range(len(f.nodes()))])\r\nplt.show()\r\n```\r\n\r\nThis message was quite brief. If something was unclear, I'd be happy to clarify those bits."}}