We numerically study the volume density probability distribution function ( n –PDF ) and the column density probability distribution function ( \Sigma –PDF ) resulting from thermally bistable turbulent flows . We analyze three-dimensional hydrodynamic models in periodic boxes of 100 pc by side , where turbulence is driven in the Fourier space at a wavenumber corresponding to 50 pc . At low densities ( n \lesssim 0.6 cm ^ { -3 } ) the n –PDF , is well described by a lognormal distribution for average local Mach number ranging from \sim 0.2 to \sim 5.5 . As a consequence of the non linear development of thermal instability ( TI ) , the logarithmic variance of the distribution for the diffuse gas increases with M faster than in the well known isothermal case . The average local Mach number for the dense gas ( n \gtrsim 7.1 cm ^ { -3 } ) goes from \sim 1.1 to \sim 16.9 and the shape of the high density zone of the n –PDF changes from a power-law at low Mach numbers to a lognormal at high M values . In the latter case the width of the distribution is smaller than in the isothermal case and grows slower with M . At high column densities the \Sigma –PDF is well described by a lognormal for all the Mach numbers we consider and , due to the presence of TI , the width of the distribution is systematically larger than in the isothermal case but follows a qualitatively similar behavior as M increases . Although a relationship between the width of the distribution and M can be found for each one of the cases mentioned above , these relations are different form those of the isothermal case .