We propose a sampling algorithm to perform system identification from a set of input-output graph signal pairs. The dynamics of the systems we study are given by a partially known adjacency matrix and a generic parametric graph filter of unknown parameters. The methodology we employ is built upon the principles of annealed Langevin diffusion. This enables us to draw samples from the posterior distribution instead of following the classical approach of point estimation using maximum likelihood. We investigate how to harness the prior information inherent in a dataset of graphs of different sizes through the utilization of graph neural networks. We demonstrate, via numerical experiments involving both real-world and synthetic networks, that integrating prior knowledge into the estimation process enhances estimation performance.