Probability density functions of the involved random variables are essential for the reliability-based design of offshore structures. The objective of present research was the derivation of probability density function (PDF) for the local joint flexibility (LJF) factor, fLJF, in two-planar tubular DK-joints commonly found in jacket-type offshore structures. A total of 162 finite element (FE) analyses were carried out on 81 FE models of DK-joints subjected to two types of axial loading. Generated FE models were validated using available experimental data, FE results, and design formulas. Based on the results of parametric FE study, a sample database was prepared for the fLJF values and density histograms were generated for respective samples based on the Freedman-Diaconis rule. Nine theoretical PDFs were fitted to the developed histograms and the maximum likelihood (ML) method was applied to evaluate the parameters of fitted PDFs. In each case, the Kolmogorov-Smirnov and chi-squared tests were used to evaluate the goodness of fit. Finally, the Inverse Gaussian model was proposed as the governing probability distribution function for the fLJF. After substituting the values of estimated parameters, two fully defined PDFs were presented for the fLJF in tubular DK-joints subjected to two types of axial loading.