Lets consider a simple derivation of GPs I've seen.
Given a dataset $ ( X, \vec{y} ) = (\vec{x}_i, y_i)_{i=0}^N $, we wish to find a function $ f $ such that $ y_i = \vec{w} \cdot \vec{\phi}(\vec{x_i}) + \epsilon_i $ where $ \vec{\epsilon} \sim \mathcal{N}\big( 0, \sigma^2 \mathbb{1} \big) $ is the gaussian noise. From various sources, it seems typical to assume $ w \sim \mathcal{N}(0, \mathbb{1}) $ (sometimes with a non-unit variance).
Then, $ \vec{y} \sim \mathcal{N}(0, \Phi \Phi^T + \sigma^2) $ with $ \Phi_{ij} = \phi_i(\vec{x}_j) $. This looks great because we don't have any parameters anymore, but this simple computes $ Pr(\vec{y} | X) $ and ignores the output of $ \vec{w} \Phi $!
I would imagine we would want to add priors for $ \phi $ and use $ Pr(\vec{y} | \Phi, X) = \mathcal{N}(w \Phi, \sigma^2 \mathbb{1}) $ (by inverting the noise distribution) for prediction.
How are we justified in ignoring $ \vec{w} \Phi $ and only using $ Pr(\vec{y}) $? What am I misunderstanding?