I’m having trouble generating out-of-sample predictions for a simple linear model with missing covariate data.

I understand how to fit a model with missing covariate data using a numpy masked array:

```
import numpy as np
import pymc3 as pm
# generate data with missing values
X_train = np.random.randn(4,2)
y_train = np.random.randn(4)
mask = [[0, 0], [1, 0], [0, 1], [0, 0],]
X_train = np.ma.masked_array(X_train, mask=mask)
# fit model
with pm.Model() as model:
alpha = pm.Normal('alpha', mu=0, sd=10)
beta = pm.Normal('beta', mu=0, sd=10, shape=2)
X_modeled = pm.Normal('X', mu=0, sd=1, observed=X_train)
mu = alpha + sum([beta[i] * X_modeled[:,i] for i in range(2)])
sd = pm.HalfCauchy('sd', beta=10)
y = pm.Normal('y', mu=mu, sd=sd, observed=y_train)
step = pm.NUTS()
trace = pm.sample(500, step=step)
```

I also understand how to generate out of sample predictions using theano shared variables when no data is missing:

```
import numpy as np
import pymc3 as pm
import theano
# generate train data
X_train = np.random.randn(4,2)
y_train = np.random.randn(4)
X_input = theano.shared(X_train)
# fit model
with pm.Model() as model:
alpha = pm.Normal('alpha', mu=0, sd=10)
beta = pm.Normal('beta', mu=0, sd=10, shape=2)
mu = alpha + sum([beta[i] * X_input[:,i] for i in range(2)])
sd = pm.HalfCauchy('sd', beta=10)
y = pm.Normal('y', mu=mu, sd=sd, observed=y_train)
step = pm.NUTS()
trace = pm.sample(500, step=step)
# generate predictions on test data with fit model
X_test = np.random.randn(4,2)
X_input.set_value(X_test)
ppc = pm.sample_ppc(trace, model=model)
```

I can’t seem to do both at the same time. Is there any way to fit a linear model with missing data and subsequently generate out of sample predictions using that fitted model with new data (which may also contain missing values)?