Categorical model with continuous dependent variable

Questions
1

I am having trouble setting up a Categorical model with the dependence that interests me.

For background, I have a probability which I observe to be an exponentially decaying function of some continuous observable parameter b. So I have a series [b, s] where I observe a binary state s given the continuous parameter b. I am not sure that a logistic regression will correctly capture the dependence on b that I am seeing, so I want to try another approach.

So I want to try modeling the probability as a beta function with a mean that depends on b.

Regardless of whether or not this will be an effective model, I should be able to set it up in pymc3, but I’m running into trouble. In reality I have thousands of samples, which is why I would rather have a continuous dependence on b rather than a bucketed dependence.

Here is a toy example of what I’m trying to do:

import numpy as np
import pandas as pd
import pymc3 as pm

b = np.array([100, 120, 150, 200, 1000, 10000, 15000])
s = np.array([1, 0, 1, 1, 0, 0, 0])
s = pd.Categorical(s).codes

model = pm.Model()
with model:
    sd = pm.HalfNormal('sd', sd=1)

    l1 = pm.HalfNormal('l1', sd=100)
    drift = tt.exp(-l1*b)

    p1 = pm.Beta.dist(balance_drift, balance_sd)

    obs = pm.Categorical('obs', p=tt.stack([p1, 1.-p1]), observed=s)

The error I observe is

TypeError: unsupported operand type(s) for -: 'float' and 'Beta'
2

Do you mean beta function or beta distribution? A log-beta function is available as pm.distributions.dist_math.betaln. If you are trying to get a beta distribution you should do p1 = pm.Beta('p', balance_drift, balance_sd)

3

A beta distribution.

When I make the change that you suggest, the error I get is

The error when converting the test value to that variable type:
Wrong number of dimensions: expected 0, got 1 with shape (7,).

This shows up when constructing ‘obs’, but seems to be coming from the drift term.
Do I need to transpose my data? But it is in an array. And I don’t want a separate Beta for every sample.

4

If s only has binary value 0 and 1, you should model it as a Bernoulli:
obs = pm.Bernoulli('obs', p=p1, observed=s)

5

I wanted to start with a simple case, so s is only binary right now. I will add more states later. With Bernoulli, I am still getting the same error though. Here is my new test case, with the changes suggested.

import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as tt

b = np.array([100, 120, 150, 200, 1000, 10000, 15000])
s = np.array([1, 0, 1, 1, 0, 0, 0])
s = pd.Categorical(s).codes

model = pm.Model()
with model:
    l1 = pm.HalfNormal('l1', sd=100)
    drift = tt.exp(-l1*b)
    sd = pm.HalfNormal('sd', sd=1)

    p1 = pm.Beta('p1', drift, sd)

    obs = pm.Bernoulli('obs', p=p1, observed=s)
6

Have you tried setting shape=7 on p1?

7

@gBokiau is right you should set the shape=7 on p1.

However, since you said earlier:

You should model p differently and had it multiple with b later on.

8

Hmmm. Is this a fundamental math problem, or just a limitation of pymc3?

So suppose I have many samples of some substance which changes state spontaneously. Like time bombs. All my samples have different ages. I want to estimate the half-life of my samples. So I have a series [a, s] where age is a continuous variable, and s is either “live” or “exploded”. How should I set up the calibration of the half-life, in pymc3?

One thought is to define a new quantity c(A) = count(s = 1 & a < A)/count(a < A), and fit that directly. ?

I also know that logistic regression is often used for modeling probabilities, but it’s not clear to me that it would be a natural model for a half-life (or, what I’m after next, a linear combination of half-lives).

9

In your setting, you are trying to find a function which map f(a) --> p, and pi(s|a) = Bernoulli(p) = Bernoulli(f(a))

In a logistic regression setting, f(a) is a logistic transformation of a linear function X*beta
In your example above, f(a) is a stochastic function that maps a exp(-l1*b) and sd into p. Since exp(-l1*b) takes a array as input, the final output p1 is also an array (1d tensor)

I hope this clear things up for you a bit.

10

Right, except it’s f(b). Sorry for renaming b -> a.

If I really do need a separate Beta for every sample, that is a lot of betas!

11

Yes that’s a lot of Betas, usually, you should infer the parameter in the function f() (l1 and sd), instead of the output p.