I want to implement a Gaussian Process model that is equivalent to Bayesian linear regression. This is my code. Assume that my data is 2D.

```
with pm.Model() as self.model:
l = pm.Gamma("l", alpha=2, beta=1)
nu = pm.HalfCauchy("nu", beta=5)
nu1 = pm.HalfCauchy("nu1", beta=5)
cov = nu1 ** 2 + nu ** 2 * pm.gp.cov.Linear(X.shape[1], l)
gp = pm.gp.Marginal(cov_func=cov)
sigma = pm.HalfCauchy("sigma", beta=1)
_y = gp.marginal_likelihood("y", X=X, y=y, noise=sigma)
map_trace = [pm.find_MAP()]
```

But when I try to plot the function estimated by above Gaussian process model, I see some non-linear behaviors as well.

Above figure is generated by using the output from the linear GP model. As you can see, it is not linear (may be it is a combination of 3 linear functions). Is it usually behavior of linear GP?

Thanks in advance…

Edit:

Below is the output using Linear regression. Which is the expected output.