I had this question at the back of my mind for a while. Consider any data for linear regression problem. Optimization algorithm calculates the coefficients of each feature and stops when cost function output reaches global minimum. The coefficients are single values for a given run of data.
When there are no multiple values for coefficients, how is a t-distribution generated and how are the standard error, t-statistic and p-value arrived at?
Assuming unbiasedness of $\hat{\beta}_{OLS}$ and invertibility of $X'X$, we can derive the variance of the OLS estimator as:
\begin{align} \sigma^{2}_{\hat{\beta}_{OLS}}&=[\hat{\beta}_{OLS}-E(\hat{\beta}_{OLS})][\hat{\beta}_{OLS}-E(\hat{\beta}_{OLS})]'\\ &=[\hat{\beta}_{OLS}-\beta][\hat{\beta}_{OLS}-\beta]'\\ &=E[(X'X)^{-1}X'U][(X'X)^{-1}X'U]\\ &=E[(X'X)^{-1}X'UU'X(X'X)^{-1}]\\ &=(X'X)^{-1}X'E[UU']X(X'X)^{-1}\\ &=(X'X)^{-1}X'\sigma_{u}^2I_N X(X'X)^{-1}\\ &=\sigma_{u}^2(X'X)^{-1}X'X(X'X)^{-1}\\ &=\sigma_{u}^2(X'X)^{-1} \end{align}
where $\beta$ is the true (population) parameter and $\sigma_{u}^2$ is the residual variance. With an estimator $\hat{\sigma}_{u}^2$, we can then obtain standard deviations $\sigma_{\hat{\beta}_{OLS}}$ from the elements of the diagonal of $\sigma_{u}^2(X'X)^{-1}$ by taking their square root.
This allows us to build t-statistics:
\begin{equation} \frac{\hat{\beta}_{OLS}-\beta_0}{\sigma_{\hat{\beta}_{OLS}}} \end{equation}
which are t-distributed with $N-K$ degrees of freedom. We don't need any more information to look at such a t distribution, check where on the support this point is, and calculate how much mass is beyond* that point (i.e. calculate the p-value).
*"beyond" meaning even farther away from the mean.
