How to deal with a variable size of the 'observed' variable

[viiv]

Hello,
In my model, the parameter I am trying to evaluate influences the size (length) of the array of the observed variable. Think of it as a parameter that makes a time-dependent function reach a goal faster or slower depending on that parameter, hence the results length vary a little.

As I was mentioning in my other question in my test with fake data, I have a parameter whose real values is 6000, and so I start with a Normal distribution at 5500, with an sd=1000. Thanks to @twiecki I have put pm.Metropolis(scaling=100) and this helps in searching the domain, to hopefully reach the 6000 value.

Nevertheless, when the Metropolis sampler starts, my observed variable Y has a size of 1186 (corresponding to the parameter’s value of 6000). Now when I try different values of this parameter, the size of the observation changes slightly, going to 1159 for instance, and then at the next trying point, 1164 for instance. At the second step I get the a Theano error about dimensions mismatch: ValueError: Input dimension mis-match. (input[3].shape[0] = 1159, input[4].shape[0] = 1164)

Here is how my code is setup, I tried to do the following trick at the end in my @as_op function proc_test(alpha):

    length_obs=len(obs)
    return(np.asarray(obs[:min(np.shape(Y)[1],length_obs)]))

and then the same trick for the likelihood:

# Priors for unknown model parameters
alpha = pm.Normal('alpha', mu=5500, sd=1000) #real is 6000
sigma = pm.HalfNormal('sigma', sd=1000)

# Expected value of outcome
test = proc_test(alpha)

# Likelihood of observations
Y_obs = pm.Normal('Y_obs', mu=test, sd=sigma, observed=Y[0,:min(np.shape(Y)[1],length_obs)])

# Inference
step = pm.Metropolis(scaling=100)
start={'alpha':5500, 'sigma':1}
trace = pm.sample(100,step=step,start=start)

but as I said it returns an error after the second step. I checked with the debugger and I feel like once length_obs is fixed in the likelihood definition (by the start value, at 1159 in my case), it is not really changed to 1164 like in a loop for instance.

Is there a way to deal with this kind of situation? I could not find anything so far.
Thanks!

[viiv]

Bumping my own thread, since it did not get a lot of views and no answers so far.
Is there anything I can add to give more details maybe?
I cannot believe I’m the only one facing this difficulty, or I must have a deep misunderstanding of how things should work, in which case I’ll be happy to look at a new approach to the problem

[aseyboldt]

What do you mean by “the size of the observation changes”? I’d have said the observation is an observation, it can’t change.
Could you maybe talk a bit about what you are modeling and post the whole code?

[viiv]

Hello,

Unfortunately I can’t post the model and/or code, otherwise I would have done it already, and my problem would probably have been resolved

As I mentioned, I have an external model of a given system (in C), which I am trying to invert, i.e. do parameter estimation on it. The parameter I am interested in at first is alpha.
To that end I make a call to this model in an as_op function proc_test(alpha). This function which takes the parameter of interest alpha, and produces an output that is an array (that is my observation). The observation is a time series, and therefore has a certain length. This length, corresponding to the duration of a phenomenon, is not fixed, but depends a little bit on alpha. This is why I have a variable-size observation.

What I am trying to do, and I am sure it can be done, it to use the Metropolis sampler to estimate the parameter alpha. The way I am doing it is maybe not how it should be done though, but I am trying.

I now realize I did not really copy-paste well from my code to the first post, so there was an error. Here is how I do things currently, and more details :

# my fake data (observation) is stored in a variable Y
# it has a fixed size (=1186 here) corresponding to a value of alpha=6000, the value I am trying to find with Metropolis

...

# function that calls my external model and returns the observation
@theano.compile.ops.as_op(itypes=[t.dscalar],otypes=[t.dvector])
def proc_test(alpha):
    
    ...

    # call of my external model
    subprocess.Popen("MYCODE.exe", shell=True).wait()

    #then I read the output from a file which fills 'obs'

    length_obs=len(obs)
    return(np.asarray(obs[:min(np.shape(Y)[1],length_obs)])) #my try so far at dealing with the length problem

# Estimation
with pm.Model() as model:
    # Priors for unknown model parameters
    alpha = pm.Normal('alpha', mu=5500, sd=1000) #real is 6000
    sigma = pm.HalfNormal('sigma', sd=1000)

    # Expected value of outcome
    test = proc_test(alpha)

    # Likelihood of observations
    Y_obs = pm.Normal('Y_obs', mu=test, sd=sigma, observed=Y[0,:min(np.shape(Y)[1],length_obs)]) #same trick as above for the length

    # Inference
    step = pm.Metropolis(scaling=100)
    start={'alpha':5500, 'sigma':1}
    trace = pm.sample(100,step=step,start=start)

Then it returns an error after the second proc_test evaluation (I have a counter). I checked with the debugger and I feel like once length_obs is fixed in the likelihood definition (by the start value, at 1159 in my case, for alpha=5500), it is not really changed afterwards to match the new length value generated by the new alpha that Metropolis tries (which is neither 5500 nor 6000, but somewhere close to that) and for instance it can happen that the length of observations is 1164. I then get ValueError: Input dimension mis-match. (input[3].shape[0] = 1159, input[4].shape[0] = 1164).
It seems that my trick doesn’t work because once length_obs is fixed in the likelihood definition, this definition is not updated like in a loop.

So I was wondering what kind of trick I could use to solve that?

Thanks in advance!

[aseyboldt]

Ah, I think I get it now, thanks. The observation is Y, the size of the prediction is changing, right?
I think the easiest solution would be to move the likelihood into the proc_test computation. Something along the lines of this:

def proc_test(alpha, sigma):
    # run the external program and write results into `predicted`
    ...
    observed_ = Y[:, min(Y.shape[1], len(predicted))]
    diff = (observed_ - predicted) / sigma
    return -0.5 * (diff * diff + np.log(2*np.pi) * sigma)

Then, you can use pm.DensityDist to get the likelihood into the model.

[viiv]

@aseyboldt

Ah, I think I get it now, thanks. The observation is Y, the size of the prediction is changing, right?

Exactly, thanks for the suggestion!

I input what you proposed in my proc_test function, and now I’m putting this in my with pm.Model() as model:

    # Expected value of outcome
    test = proc_test(alpha, sigma)

    # Likelihood (sampling distribution) of observations
    Y_obs = pm.DensityDist('Y_obs', logp=test, observed=Y)

But I get TypeError: 'TensorVariable' object is not callable.
I tried different syntaxes such as directly logp=proc_test or logp=proc_test(alpha,sigma) but I always get the same error. Any ideas? It feels weird to put observed=Y since we’re now looking at a sub-sample of Y in some cases for instance, but I’m not sure.
I have to admit that I am not familiar yet with DensityDist also.

[junpenglao]

You might find this post helpful where I was trying to explain the mechanism of DensityDist.

Problem model definition PyMC3

Hi, I have a custom likelihood and I’m trying to determine the posterior probabilty for my parameters , given Uniform priors. I’ve written (trying to follow https://arxiv.org/pdf/1507.08050.pdf): with Model() as paper: #priors m = Uniform('m',0.,1.) q = Uniform('q',0.,5.) sigmay = Uniform('sigmay',0.,1.) # deterministic rule logEpeak_i=(q)+(m*logEiso) #likelihood definition def likelihoodMCMC(m,q,sigmay,logEiso,logEpeak_i,logErrEiso,logErrEpeak_i): return np.exp(0.5*n…

[viiv]

Thanks,

I used the dictionary syntax for observed, and I have this slightly modified code now :

@theano.compile.ops.as_op(itypes=[t.dscalar,t.dscalar],otypes=[t.dvector])
def proc_test(Y,alpha,sigma):
    
    ...
    
    predicted= ... #what I had defined as obs before
    observed_ = Y[:, min(Y.shape[1], len(predicted))]
    diff = (observed_ - predicted) / sigma
    return -0.5 * (diff * diff + np.log(2*np.pi) * sigma)
  
with pm.Model() as model:

    # Priors for unknown model parameters
    alpha = pm.Normal('alpha', mu=5500, sd=1000) #real is 6000
    sigma = pm.HalfNormal('sigma', sd=1000)

    # Likelihood (sampling distribution) of observations
    Y_obs = pm.DensityDist('Y_obs', proc_test, observed={'Y':Y,'alpha':alpha,'sigma':sigma})

    # Inference!
    step = pm.Metropolis(scaling=100)

    start={'alpha':5500, 'sigma':1}

    trace = pm.sample(100,step=step,start=start)

But I get now a TypeError: make_node() got an unexpected keyword argument 'Y'. Not sure what is wrong again here, because it looks like I have the same syntax at all steps as the code in the thread you linked.

[junpenglao]

@viiv

@theano.compile.ops.as_op(itypes=[t.dscalar,t.dscalar],otypes=[t.dvector])

You need to specify the input for Y as well:

@theano.compile.ops.as_op(itypes=[t.dvector,t.dscalar,t.dscalar],otypes=[t.dvector])

Also, I think you should return a scaler instead of a vector:

    return np.sum(-0.5 * (diff * diff + np.log(2*np.pi) * sigma))

[viiv]

Hello,

Yes, clearly I needed to change that,
I now have @theano.compile.ops.as_op(itypes=[t.dvector,t.dscalar,t.dscalar],otypes=[t.dscalar])
nevertheless I still get the exact same error (which is why I oversaw this simple mistake)…

I did some tests, and it seems that the keyword argument in the error is the first one that is passed in observed. If I change the order of the dictionnary, by putting alpha first, it yields
TypeError: make_node() got an unexpected keyword argument 'alpha'

Going into the code of op.py it seems that the kwargs are maybe confusing make_node() here.
I tried to remove the as_op of my proc_test function, then it runs without this error, but instead of using the value of alpha it sees it as a string in the function, which of course produces errors later on.

I’m not sure if I’m using this DensityDist the right way…

[viiv]

Hello,

I looked at this make_node() error a little more, but could not find anything that would make it disappear.
I was wondering if anyone had more ideas?

Here is the code with the last remarks from @junpenglao

@theano.compile.ops.as_op(itypes=[t.dscalar,t.dscalar,t.dvector],otypes=[t.dscalar])
def proc_test(alpha,sigma,Y):
    
    ...
    
    predicted= ... #what I had defined as obs before
    observed_ = Y[:, min(Y.shape[1], len(predicted))]
    diff = (observed_ - predicted) / sigma
    return np.sum(-0.5 * (diff * diff + np.log(2*np.pi) * sigma))
  
with pm.Model() as model:

    # Priors for unknown model parameters
    alpha = pm.Normal('alpha', mu=5500, sd=1000) #real is 6000
    sigma = pm.HalfNormal('sigma', sd=1000)

    # Likelihood (sampling distribution) of observations
    Y_obs = pm.DensityDist('Y_obs', proc_test, observed={'alpha':alpha,'Y':Y,'sigma':sigma})

    # Inference!
    step = pm.Metropolis(scaling=100)

    start={'alpha':5500, 'sigma':1}

    trace = pm.sample(100,step=step,start=start)

I still get TypeError: make_node() got an unexpected keyword argument 'alpha'.
And to make it clear, if I change the order of the terms in the DensityDist dictionnary parameters, putting for instance Y in first poistion, I’ll get TypeError: make_node() got an unexpected keyword argument 'Y'.
I’m really unsure what is going on here, but it must be some kind of syntax error.

[aseyboldt]

Functions decorated with as_op can’t take keyword arguments. You could wrap it in a simple python function:

@theano...
def proc_test(alpha, sigma, Y):
    ...

with model:
    alpha = ...
    sigma = ...

    def logp(Y):
        return proc_test(alpha, sigma, Y)

    Y_obs = pm.DensityDist('Y_obs', logp, observed=Y)

[viiv]

Thanks it worked!

Here are a few more changes that needed to be done too, if anyone is interested:

observed_ = Y[:min(Y.shape[0], len(predicted))]

since I noticed it was outputting one element instead of the array. I actually noticed my Y had an extra pair of , so I removed it, and then changed the code above which returns the needed array, and then

diff = (observed_ - predicted[:min(Y.shape[0], len(predicted))]) / sigma

to match the case where I have more predictions than in Y.

Afterwards, I came across some really weird NumPy errors at about 17% in the run: TypeError:

Cannot cast ufunc multiply output from (‘float64’) to dtype(‘int32’) with casting rule ‘same_kind’

so I had to change the scale function in metropolis.py this way (it didn’t like the *= apparently), but now it works:

# Switch statement
if acc_rate < 0.001:
    # reduce by 90 percent
    np.multiply(scale, 0.1, out=scale, casting='unsafe')
elif acc_rate < 0.05:
    # reduce by 50 percent
    np.multiply(scale, 0.5, out=scale, casting='unsafe')
elif acc_rate < 0.2:
    # reduce by ten percent
    np.multiply(scale, 0.9, out=scale, casting='unsafe')
elif acc_rate > 0.95:
    # increase by factor of ten
    np.multiply(scale, 10.0, out=scale, casting='unsafe')
elif acc_rate > 0.75:
    # increase by double
    np.multiply(scale, 2.0, out=scale, casting='unsafe')
elif acc_rate > 0.5:
    # increase by ten percent
    np.multiply(scale, 1.1, out=scale, casting='unsafe')

By the way, one last thing @aseyboldt , it seems to me that in your reply where you tell me to return the logp in the function directly, you write to put

return -0.5 * (diff * diff + np.log(2*np.pi) * sigma)

which @junpenglao changed into

return np.sum(-0.5 * (diff * diff + np.log(2*np.pi) * sigma))

but shouldn’t it be

return np.array(np.sum(-0.5 * (diff * diff + np.log(2np.pi) + 2np.log(sigma))))

instead? (nevermind the np.array that was also creating some bugs).
If you do the math for the log of the Normal pdf:

log[1 / (sigma sqrt(2pi)) * exp(-0.5 * (x-mu)^2 / sigma^2)] = -0.5 * [(x-mu)^2 / sigma^2 + log(2pi) +2*log(sigma)]

so the last term is not * sigma, but +2*log(sigma), no?

[aseyboldt]

You are right about the logp of the normal, sorry about that. One would think that I’ve used that enough times to get it right…

About the dtype issue: Could you try to change the scaling to 100.0 instead of 100? I think this a bug in pymc, where we don’t check that the scaling in a floating point number.

[viiv]

I will try that, revert my changes in metropolis.py, and report if there are any bugs, thanks for the suggestion!

[viiv]

I can now confirm that I had no problem with that during the whole run.