The paper focuses on the numerical approximation of nabla fractional order systems with the conditions of nonzero initial instant and nonzero initial state. First, the inverse nabla Laplace transform is developed and the equivalent infinite dimensional frequency distributed models of discrete fractional order system are introduced. Then, resorting the nabla Laplace transform, the rationality of the finite dimensional frequency distributed model approaching the infinite one is illuminated. Based on this, an original algorithm to estimate the parameters of the approximate model is proposed with the help of vector fitting method. Additionally, the applicable object is extended from a sum operator to a general system. Three numerical examples are performed to illustrate the applicability and flexibility of the introduced methodology.