The temperature dependencies of NMR shifts in the critical region of two coexisting phases have been simulated using statistical thermodynamics and graph-theory consideration of equilibrium processes of molecular association. Microparameters of magnetic screening of various water and water/pyridine structures used in the statistical averaging have been evaluated by density functional theory calculations (PBE1PBE and B3PW91 functionals in the 6-311++G** basis set). The gauge-including atomic orbital (GIAO) approach has been applied to ensure gauge invariance of the results. Solvent effects were taken into account by a polarized continuum model (PCM). NMR shifts "order parameters" (Delta delta = parallel to delta(+) - delta(-)parallel to) and "diameters" (phi delta = parallel to(delta(+) + delta(-))/2 -delta(C)parallel to, where delta(+), delta(-), and delta(C) are the chemical shifts of coexisting phases and at the critical point respectively) have been calculated in each case close to the lower critical solution point (T-L) and processed using linear regression analysis of Delta delta similar to parallel to T - T-L parallel to and phi delta similar to parallel to T - T-L parallel to in the log-log plot. It has been shown that the critical index beta can be evaluated with high precision from the slope of Delta delta = f(T - T-L) at any realistic set of model input parameters. The slope of diameter has been found to depend on both input beta and alpha values. The obtained phi delta slopes (0.58-0.63) are very close to 2 beta values. The results are discussed within the concept of complete scaling. Results of simulation are compared and supported by experimental NMR data for water/2,6-lutidine, acetic anhydride/n-heptane, and acetic anhydride/cyclohexane systems.