Predicting on hold-out data，Error：Input dimension mis-match. (input[1].shape[0] = 126, input[2].shape[0] = 3)

Questions
1

In a simple model with 42 data, I can predict on hold-out data. but In another complex model with 126 data, the model throw the error:
’’ ValueError: Input dimension mis-match. (input[1].shape[0] = 126, input[2].shape[0] = 3)
Apply node that caused the error: Elemwise{Composite{exp((i0 + (i1 * i2) + (i3 * i4) + i5))}}[(0, 1)](InplaceDimShuffle{x}.0, AdvancedSubtensor1.0, <TensorType(int64, vector)>, InplaceDimShuffle{x}.0, <TensorType(int64, vector)>, AdvancedSubtensor1.0)
Toposort index: 5
Inputs types: [TensorType(float64, (True,)), TensorType(float64, vector), TensorType(int64, vector), TensorType(float64, (True,)), TensorType(int64, vector), TensorType(float64, vector)]
Inputs shapes: [(1,), (126,), (3,), (1,), (3,), (126,)]
Inputs strides: [(8,), (8,), (8,), (8,), (8,), (8,)]
Inputs values: [array([-0.]), ‘not shown’, array([5, 6, 7], dtype=int64), array([-0.]), array([40, 40, 40], dtype=int64), ‘not shown’]
Outputs clients: [[‘output’]]

HINT: Re-running with most Theano optimization disabled could give you a back-trace of when this node was created. This can be done with by setting the Theano flag ‘optimizer=fast_compile’. If that does not work, Theano optimizations can be disabled with ‘optimizer=None’.
HINT: Use the Theano flag ‘exception_verbosity=high’ for a debugprint and storage map footprint of this apply node. ‘’

My model like this:

x_shared = shared(elec_year)
x_shared1 = shared(elec_tem)
y_shared = shared(elec_faults)

with pm.Model() as unpooled_model:

    sigma = pm.HalfCauchy('sigma', 8, testval=1.)
    beta = pm.Normal('beta', 0, 10)
    beta1 = pm.Gamma('beta1', 2, 10, shape=companiesABC)
    beta2 = pm.Normal('beta2', 0, 10)
    sita = pm.Gamma('sita', 1, 100000)
    u = pm.Normal('u', 0, sita, shape=companiesABC)

    mu = pm.Deterministic('mu', tt.exp(beta + beta1[companyABC]*x_shared + beta2*x_shared1 + u[companyABC]))

    Observed = pm.Weibull("Observed", alpha=mu, beta=sigma, observed=y_shared)  
    start = pm.find_MAP()
    trace2 = pm.sample(4000, start=start)

x_shared.set_value([5, 6, 7])
x_shared1.set_value([40, 40, 40])
y_shared.set_value([0, 0, 0])

with unpooled_model:
    post_pred = pm.sample_ppc(trace2, samples=500)
post_pred['Observed'].mean(axis=0)

whether I have to change the x_shared.set_value([5, 6, 7]) into a 126 length data?

1 Like
2

The code syntax seems fine to me, I would suggest you to double check the data input.

3

Thanks, I tried the 126 length data, and it works well

4

you have answered this question:grinning:
however, I run this example and my another question is why the predicted value is not correct in this example.

5

A few points for you to further consider, assumpting by not correct you meant the 95% prediction does not include the actual observation:
1, check the model fitting, is it converged? any divergence? is the trace look fine? what about the rhat and effective sample size.
2, if possible, compare with a frequentistic solution. since in this case it is a GLM you can compare it with the solutions from statsmodels or numpy.linalg.solve
3, if the model fitting is fine, and the estimation is similar to a frequentistic solution, then it would means the model is not good enough. Then try implementing better model (by accounting for covariance), use other more robust option (e.g., student t intead of Gaussian to account for outliners) etc.

6

Thanks so much much again, I will try it:grinning:

7

Hi, elder brother peng, do you know how to use python achieve the maximum likelihood estimation ? I want to do some comparison but I don’t know whether I can put MLE into model of Pymc3(may can not) .If not ,than I have to make another model in order to use MLE? I have so many quesitons… :joy:

8

LOL, do you mean doing MLE using say scipy on a linear model, and then plug the estimation into a PyMC3 model to compare prediction? You surely can. For the first part you can do a search online there are quite a lot of tutorial. As for plugging value into a PyMC3 model and do prediction, you can follow a workflow similar to below:

Say you generated some data following one of the GLM doc (http://docs.pymc.io/notebooks/GLM-linear.html). Now you build your model using the glm module:

with Model() as model:
    glm.GLM.from_formula('y ~ x', data)

This model has logp function that takes a dict as input and output loglikelihood, the input format is:

model.test_point
Out[4]: {'Intercept': array(0.0), 'sd_log__': array(0.0), 'x': array(0.0)}

Notice here that sd is in log space, here is an explanation http://docs.pymc.io/notebooks/api_quickstart.html#Automatic-transforms-of-bounded-RVs. And x is the slope.

Now you can feed a dict with a similar structure as above to evaluate logp, for example:

model.logp(model.test_point)

or even feed the real value:

model.logp(dict(Intercept=true_intercept,sd_log__=np.log(.5),x=true_slope))

And you can find the MLE of model:

with model:
    map1 = find_MAP()
    
logp = -170.67, ||grad|| = 44.491: 100%|██████████| 21/21 [00:00<00:00, 544.33it/s]  

In[15]: map1
Out[15]: 
{'Intercept': array(0.9488320028422407),
 'sd': array(0.48009546387935403),
 'sd_log__': array(-0.7337703117728311),
 'x': array(2.0195327989981764)}

For PPC, you can past a list containing the MLE to sample_ppc:

ppc = sample_ppc([map1], samples=100, model=model)

I think it might be helpful to follow a book to learn Bayesian Statistics with PyMC3. You can find PyMC3 ports of two great Bayesian book here: https://github.com/pymc-devs/resources, which I encourage you to follow.

9

I do not know what to say but thanks… :smile: I am studying the book that you program for it . :star_struck:

1 Like