Is the wrong answer just because of the mean alignment?
Admittedly, I have only worked with your cherry blossom example, but when I input new year values extended out and just far enough back in time to maintain the mean year… I get (by the eyeball) very similar results.
from your git hub (using: phi, sqrt_psd = gp.prior_linearized(X) )
After alteration:
It could very well be quirk that this happened to work out… just curious about what would lead to your intuition about getting the wrong answer from a projection using pm.set_data on a HSGP.
What I did instead:
x = df[“year”].values
with pm.Model() as model:
model.add_coord(‘year’, df.year.tolist(),mutable=True)
eta = pm.Exponential(“eta”, lam=0.25, shape=1)
ell_params = pm.find_constrained_prior(
pm.InverseGamma,
lower=20, upper=200,
init_guess={“alpha”: 2, “beta”: 20},
mass=0.95,
)
ell = pm.InverseGamma(“ell”, **ell_params)
# Specify covariance function and GP
cov = eta[0]**2 * pm.gp.cov.Matern52(1, ls=ell)
intercept = pm.Normal(“intercept”, mu=df[“temp”].mean(), sigma=2 * df[“temp”].std())
gp = pm.gp.HSGP(m=[400], L=[1565.0], cov_func=cov, mean_func=pm.gp.mean.Constant(intercept))
x_mu = pm.ConstantData(‘x_mu’,np.mean(x))
X = pm.MutableData(f"X_mu_{x_mu}", x[:, None])
X = X
f = gp.prior(“f”, X=X, dims=“year”)
mu = pm.Deterministic(“mu”, f, dims=“year”)
sigma = pm.HalfNormal(“sigma”, sigma=2 * df[“temp”].std())
pm.Normal(“y”, mu= f, sigma=sigma, observed=df[“temp”].values, shape=X.shape[0 ])
Then:
def extend_x_keep_mean(x_extension=50, target_mean=np.mean(x), x=x, make_int=True):
maintain_mean_of_x_halfwidth = x_extension + np.max(x) - target_mean
x_new = np.arange((target_mean - maintain_mean_of_x_halfwidth),
(target_mean + maintain_mean_of_x_halfwidth))
if make_int:
x_new = x_new.astype(int)
print(f’{x_new.mean()-target_mean=}')
return x_new
x_new = extend_x_keep_mean()
Then:
with model:
model.set_dim(“year”, new_length=len(x_new),coord_values=x_new.tolist())
model.set_data(name=f"X_mu_{x_mu}",values=x_new[:, None], coords={“year”: x_new.tolist()})
ppc = pm.sample_posterior_predictive(idata_gp, var_names=[ ‘f’])