Composing partial models?

[burnpanck]

When experimenting with complex hierarchical models, I try to build the model in a modular fashion, with a clean structured interface between the weakly coupled parts of the model. While not much of an actual problem in practice, I find it conceptually unfitting that the modular parts need to provide a reference to the complete Model while it is constructed. In my view, a model is just a collection of random variables that are either dependent on other random variables, or free random variables that are associated with a probability distribution. Such collections should be easy to compose. Thus I'm wondering if there is an API to compose multiple related or unrelated models.

In fact, I think it should be possible to specify a model by simply describing the relations between random variables without ever explicitly referencing a global model with which the variables are associated. In that view, the model could in principle always stay implicit, as a user you simply ask PyMC3 to calculate/sample a posterior distribution of a given set of random variables, and PyMC3 would then traverse the theano graph to identify all involved random variables.

An intermediate version of these two extremes would be an API that allows a user to explicitly ask PyMC3 to generate a model by collecting all involved random variables from a set of random variables of explicit interest, thus removing the reference to the model during construction of the random variables. Indeed, as far as I can tell from reading the code, the random variables themselves don't do much with the model that they are so eager to know at construction time, except storing it in an attribute. The model on the other hand collects the variables directly, but it shouldn't be too difficult to delay this collection effort to later.

I am inclined to try to implement a mock Model which doesn't store anything to keep the random variables code happy, and then add a function to collect random variables into a model at a later time. Do you see any difficulties with that approach, i.e. are there situations where the random variables actually need to know their model at construction time?

[twiecki]

Related: #1525 CC @ferrine.

[twiecki]

In general, that sounds useful. Can you first write some pseudo-code of what your dream API might look like and what it would allow?

[burnpanck]

I described three API's among which the third used to be my favourite (the compromise).
However, while writing, I realised that I possibly need to think a bit more about the later two.
Here is how I imagined them from the user's point of view:

1. Composing existing models:

with Model() as A:
    a = Normal('a',1,0)
    adata = Normal('adata',a,1,observed=-1)
with Model() as B:  # completely unobserved
    b = Normal('a',1,0)
with Model() as Parent:
    parent_data = Normal('parent_data',a*b,1,observed=2)
full_model = Model.combine_models(A,B,Parent)
trace = sample(100000,model=full_model)

Of course, the example is silly, but just to illustrate the proposed API, which is the second-to-last line:
A simple function to generate a model that is the union of all random variables of the submodels.
However, there are a few things to note: I don't know if the current code would accept B, as there is no data. Also, when defining Parent, I refer to random variables from "other" models, which might lead to problems with the current implementation.

2. Completely implicit models:

a = Normal('a',1,0)
amsmt = Normal('amsmt',a,1)
b = Normal('a',1,0)
parent_msmt = Normal('parent_msmt',a*b,1)
trace = sample(100000,infer=[a],data={amsmt:adata,parent_msmt:parent_data})

Here, no model was ever explicitly specified, but when PyMC3 gets asked to infer the value of a, it should implicitly collect all random variables that are related to a. Which, as I write this, raises an important point: From the point of view of the theano graph, a is not a function of any of the data, therefore we need to explicitly associated the actual data with the model. Of course, we could also have used observed= to associated the data, but we would still have to explicitly mention the observed random variables, as they do not feature in the ancestry graph of a.
Another problem here is that we will have to do these specifications again whenever we want to actually work with the model/data, (i.e. find_MAP or whatever). Therefore the compromise:

3. Compromise

a = Normal('a',1,0)
amsmt = Normal('amsmt',a,1)
b = Normal('a',1,0)
parent_msmt = Normal('parent_msmt',a*b,1)
model = Model.collect_RVs(vars=[a],data={amsmt:adata,parent_msmt:parent_data})
trace = sample(100000,model=model)

Same issues as with API 2.

Which leaves only API 1 as a simple option for me now.
I should probably have a look at #1525, and think some more.

[fonnesbeck]

The important constraint is that we never want the user to have to explicitly add variables to a model, as you see in several of the ML libraries (i.e. model.add(variable)). In PyMC2, you could declare a set of variables that were automatically discovered in the current environment by the Model object, which was a bit of a hack, but convenient. Personally, I like the use of context managers, since they are an established Python idiom, yet they provide enough "magic" to automate model building to the degree that I like. However, I share your enthusiasm for being able to combine models, or create arbitrary sub-models; I am a big fan of multi-model inference, and being able to easily define a set of nested models appeals to me.

[ferrine]

I'm currently working on good api for convenient usage of different models. Now I'm not able to do:

with Model() as A:
    a = Normal('a',1,0)
    adata = Normal('adata',a,1,observed=-1)
with Model() as B:  # completely unobserved
    b = Normal('a',1,0)
with Model() as Parent:
    parent_data = Normal('parent_data',a*b,1,observed=2)
full_model = Model.combine_models(A,B,Parent)
trace = sample(100000,model=full_model)

But I can

full_model = Model()
# choose name for model or else names will overlap 
with Model(name='suba', model=full_model) as A:
    a = Normal('a',1,0)    # name will be 'suba_a'
    adata = Normal('adata',a,1,observed=-1)
with Model(name='subb',model=full_model) as B:  # completely unobserved
    b = Normal('a',1,0)    # name will be 'subb_a'
with Model(model=full_model) as Parent:
    parent_data = Normal('parent_data',a*b,1,observed=2)
trace = sample(100000,model=full_model)

Note: it was possible to run with Model(name='suba'): ... in the context of full_model to get the same effect
Is it simple api or you have thoughts how I can make it better?

[burnpanck]

@ferrine: My goal was not having to declare full_model beforehand, because model suba really doesn't care at all about full_model so why the need to mention it? Of course, when the code is as linear as in this simple example, we could as well just have wrapped everything in a single model.

@fonnesbeck: I completely agree, the user should not have to add variables to the model. I certainly do not want to imply that the context manager should go away. However, I also think that it should not be necessary to declare the model beforehand, to simplify separation of concepts when building complex multi-models. I was trying to come up with an API that worked without the context manager, because I wanted to point out that it should be possible to build it simply from variables that could have been created beforehand. But since at some point, one will need to group the variables, so why not do so right when declaring them, bringing me back to API 1.

Maybe one comment to the context manager syntax: One other option would have been to instead use a class declaration. This would fit better my view of a model as a general description of relationships between random variables, and instances of that model could be a particular association with data or (random) variables external to the model. With enough meta-class black magic, one would certainly be able to achieve most of the features found on the context manager. One advantage would be that the names of the variables would not have to be repeated. A disadvantage would be that variables declared deeply inside nested scopes would not automatically be captured by the model. I particularly see merit of such an approach again with "lasagne" models, where nesting would then rather occur in the form of actual sub-models that get passed on to the main model, such that all related variables are properly found by the main model.

[ferrine]

@burnpanck
see PR #1525 for class based models tutorial and sample implementation of LinearComponent and Glm. This will be merged soon I hope. Briefly it looks like this:

class CustomModel(Model):
    # 1) override init
    def __init__(self, mean=0, sd=1, name='', model=None):
        # 2) call super's init first, passing model and name to it
        # name will be prefix for all variables here
        # if no name specified for model there will be no prefix
        super(CustomModel, self).__init__(name, model)
        # now you are in the context of instance,
        # `modelcontext` will return self
        # you can define variables in several ways
        # note, that all variables will get model's name prefix

        # 3) you can create variables with Var method
        self.Var('v1', Normal.dist(mu=mean, sd=sd))
        # this will create variable named like '{prefix_}v1'
        # and assign attribute 'v1' to instance
        # created variable can be accessed with self.v1 or self['v1']

        # 4) this syntax will also work as we are in the context
        # of instance itself, names are given as usual
        Normal('v2', mu=mean, sd=sd)

        # something more complex is allowed too
        Normal('v3', mu=mean, sd=HalfCauchy('sd', beta=10, testval=1.))

        # Deterministic variables can be used in usual way
        Deterministic('v3_sq', self.v3 ** 2)
        # Potentials too
        Potential('p1', tt.constant(1))

# After defining a class CustomModel you can use it in several ways

# I:
#   state the model within a context
with Model() as model:
    CustomModel()
    # arbitrary actions

# II:
#   use new class as entering point in context
with CustomModel() as model:
    Normal('new_normal_var', mu=1, sd=0)

# III:
#   just get model instance with all that was defined in it
model = CustomModel()

# IV:
#   use many custom models within one context
with Model() as model:
    CustomModel(mean=1, name='first')
    CustomModel(mean=2, name='second')

[burnpanck]

Indeed, I already suspected that you have some sweetness in #1525, though I didn't have time yet to understand it fully. If time permits, I will try to reformulate my model in terms of this patch, and see how it fits my concept. In any case, thanks to all of you for the discussion and of course for PyMC3!

[twiecki]

@burnpanck Yes, would be great if you test #1525 (it's on the 3.1 branch now).

[twiecki]

@burnpanck As we merged #1524, can we close this or are some things not supported?

[burnpanck]

I have to admit, I did not get around to test it yet. I had succeeded in modelling my problem without modular composition, and since I'm now very busy for the next two months, it is very unlikely that I will play with this any time soon. So I think closing is fine. Thanks again for your support!