# Tag Info

15

You can compute the likelihood of a model indeed using model.logp(). As input, it requires a point. For example, the BEST model from the examples directory I can do: np.exp(model.logp({'group1_mean': 0.1, 'group2_mean': 0.2, 'group1_std_interval': 1., 'group2_std_interval': 1.2, ...

14

It's generally true in my personal experience as a professional data scientist. It's true in my personal experience because it's what I observe most of the time. If you're asking why it happens this way, it's for a few reasons: Many traditional ML algorithms are nowadays available "off the shelf", including sophisticated ensemble methods, neural networks, ...

13

You are looking the what's called the predictive distribution. To include this is very simple. Before creating the Model, add the additional stochastic variable: predictive = mc.Normal( "predictive", mean, precision ) model = Model( {"mean": mean, "obs": obs, "pred":predictive}) ... predictive_traces = mcmc.trace("predictive")[:] hist( predictive_traces ) ...

12

This is of course a diverse set of people with a range of opinions getting together and writing a wiki. I summarize I know/understand with some commentary: Choosing your prior based on computational convenience is an insufficient rationale. E.g. using a Beta(1/2, 1/2) solely because it allows conjugate updating is not a good idea. Of course, once you ...

11

Are you absolutely certain that half came from one distribution and the other half from the other? If not, we can model the proportion as a random variable (which is a very bayesian thing to do). The following is what I would do, some tips are embedded. from pymc import * size = 10 p = Uniform( "p", 0 , 1) #this is the fraction that come from mean1 vs ...

9

A similar problem is discussed in Gelman, Bayesian Data Analysis, (2nd ed, p. 128; 3rd edition p. 110). Gelman suggests a prior $p(a,b)\propto (a+b)^{-5/2}$, which effectively constrains the "prior sample size" $a+b$, and therefore the beta hyperprior is not likely to be highly informative on its own. (As the quantity $a+b$ grows, the variance of the beta ...

9

Assuming your likelihood and MCMC work correctly, which we can't know of course, and assuming your new variables are not ineffective: If you increase the number of parameters under calibration with data fixed, an increase in marginal posterior uncertainty up to the point that the marginal posterior looks like the prior is quite common. This happens if your ...

8

When defining w, the p parameter must be a list of doubles, not a list of lists of doubles. This means you have to define a w variable for each word in each document. Also it helps to 'complete' the Dirichlet variables using the CompletedDirichlet function. Here is the working code: import numpy as np import pymc as pm K = 2 # number of topics V = 4 # ...

8

Please use with PyMC 2.3.2 It looks like there is a bug in 2.3.4 (the most recent version) that is causing the wrong inference. This took me a while to discover, but it was solved when I downgraded to PyMC 2.3.2. Model: import pymc as pm p = [ #brown, yellow, red, green, orange, tan, blue [.3, .2, .2, .1, .1, .1, .0 ], # 1994 bag [.13, .14, .13, ....

8

They do not provide any scientific/mathematical justification for doing so. Most of the developers do not work on this kind of priors, and they prefer to use more pragmatic/heuristic priors, such as normal priors with large variances (which may be informative in some cases). However, it is a bit strange that they are happy to use PC priors, which are based ...

7

I would use a latent variable approach, since that's what x an y are. However, its not clear that all four parameters would be identifiable in this case. It would be helpful if you had some prior information for one or two of them. Here's an example: import pymc as pm # Priors mu_x = pm.Normal('mu_x', 0, 0.001, value=0) sigma_x = pm.Uniform('sigma_x', 0, ...

7

This does seem like a good model, implemented correctly in PyMC. There are two Bayesian stats facts that we can use to confirm your answer with another method: $\textrm{Beta}(1,1)$ is equivalent to the uniform distribution on the interval $[0,1]$; The beta and binomial distributions are conjugate. This means that the posterior distribution of $p$ is also ...

7

I'll record my modelling thoughts as I create it: import pymc as pm How can I recreate this data? Well, N_T cars enter. Do I know N_T? No, so it's a random variable. For simplicity, I'll say its a discrete uniform with max 1000 (you can change this to, say, a Poisson) N_T = pm.DiscreteUniform('N_T', 0, 1000) Of these N_T cars that enter, I know a car ...

7

The way I understand it, PGMs is a broad class including Markov Random Fields, Conditional Random Fields and Bayesian Networks (to name a few). PyMC works on Bayesian Networks (i.e. if the network can be represented as a Directed Acyclic Graph).

7

Author here. There are a few points I can elaborate on: Your data is going to be distributed across many computers in a real cluster, so each computer has a fraction of the data. Locally, then, the computations should run fine. But the question next is: how do you combine all these inferences? Sure you could just merge the traces from each computer, but ...

7

We use pm.Potential here primarily to get around the definition of a likelihood. We ordinarily use it to constrain our likelihood in the manner described in the PyMC docs, but in this example we never end up defining a true likelihood (which would require the inclusion of observations). As such, all the samples that we draw are based on how we defined the ...

6

Take a look at a post in Healthy Algorithm: http://healthyalgorithms.com/2011/11/23/causal-modeling-in-python-bayesian-networks-in-pymc/ also in PyMC's totorial: http://pymc-devs.github.io/pymc/tutorial.html Maybe you would try the following code clip (assuming you have imported pymc as mc): A = mc.Normal('A', mu_A, tau_A) B = mc.Normal('B', mu_B, tau_B) ...

6

A couple of points, related to the discussion above: The choice of diffuse normal vs. uniform is pretty academic unless (a) you are worried about conjugacy, in which case you would use the normal or (b) there is some reasonable chance that the true value could be outside the endpoints of the uniform. With PyMC, there is no reason to worry about conjugacy, ...

6

SeanEaster has some good advice. Bayes factor can be difficult to compute, but there are some good blog posts specifically for Bayes factor in PyMC2. A closly related question is goodness-of-fit of a model. A fair method for this is just inspection - posteriors can give us evidence of goodness-of-fit. Like quoted: "Had no change occurred, or had the ...

6

I hope you like Python! I'll recite my comment here: This sounds like a hierarchical model. If I wanted to recreate the dataset, here's what I'd do: Let $D$ be a $Beta(\alpha, \beta)$ distribution (reasonable since we are dealing with probabilities). We don't know $\alpha, \beta$, we assign priors to them, say exponential for both with some $\lambda$ ...

6

The problem is caused by the way that PyMC draws samples for this model. As explained in section 5.8.1 of the PyMC documentation, all elements of an array variable are updated together. For small arrays like center this is not a problem, but for a large array like category it leads to a low acceptance rate. You can see the acceptance rate via print mcmc....

6

From a pure implementation perspective, it should be straightforward: take your model code, replace every trainable Variable creation with ed.Normal(...) or sth similar, establish variational posteriors as well, zip them in a dict, feed it to some inference object from edward et voila. The problem is that variational training of RNNs, since based on ...

5

Gelman has good advice for setting priors for variance parameters in Bayesian models. There is too much structure in this model for the data you are trying to fit. In particular, it is not clear why you are modeling mu, rather than just putting a prior on it. The way you have it set up, you are claiming that mu is sampled from another normal with unknown ...

5

You could try the approach recommended by Steve Goodman and calculate the minimum bayes factor: Toward Evidence Based Medical Statistics 2: The Bayes Factor To get this from mcmc results, you can subtract the estimate for the group level parameters for each step to get a posterior distribution of the difference as was done by John Kruschke in this paper: ...

5

I think there are a few approaches here. First Approach As far as I know, there is no way to use @deterministic or @stochastic (without the likelihood). An alternative way is to use the potentials class, which is like multiplying your likelihood by a factor. In this case, we should multiply by the pdf of a lognormal given $Z$ and $X$. import pymc as mc ...

5

This took quite a bit of work, but I got it in the end. Note that I used the development version (pymc 2.2grad) from github, not the older version available on pypi. Also, this runs pretty slowly, since it doesn't make good use of numpy's array manipulation. For data sets of reasonable size, some smart preprocessing and a rewrite of the cutpoints node ...

5

PyMC2 can be combined with the LP solver of your choice to solve stochastic LP problems like this one. Here is code to do it for this very simple case. I've left a note on how I would change this for a more complex LP. c1 = pm.Normal('c1', mu=2, tau=.5**-2) c2 = -3 b1 = pm.Normal('b1', mu=0, tau=3.**-2) @pm.deterministic def x(c1=c1, c2=c2, b1=b1): # ...

5

That's an interesting problem. I think it nicely highlights the power of Bayesian statistics as you can build these custom models based on how you think your data was generated and then invert them. To answer your question: I would suggest to instead do the mixture inside of the likelihood where you can choose the specific likelihood function for each ...

5

Autocorrelation dictates the amount of time you have to wait for convergence. If autocorrelation is high, you will have to use a longer burn-in, and you will have to draw more samples after the burn-in to get a good estimate of the posterior distribution. Low autocorrelation means good exploration. Exploration and convergence are essentially the same ...

5

The term "hierarchical" in this example means that it is a hierarchical Bayesian model. It is not a hierarchical GLM in the sense you describe.

Only top voted, non community-wiki answers of a minimum length are eligible