We provide a framework to estimate the systematic uncertainties in chemical freeze-out parameters extracted from $chi^2$ analysis of thermal model, using hadron multiplicity ratios in relativistic heavy-ion collision experiments. Using a well known technique of graph theory, we construct all possible sets of independent ratios from available hadron yields and perform $chi^2$ minimization on each set. We show that even for ten hadron yields, one obtains a large number ($10^8$) of independent sets which results in a distribution of extracted freeze-out parameters. We analyze these distributions and compare our results for chemical freeze-out parameters and associated systematic uncertainties with previous results available in the literature.