Derivation of Box-Cox and Yeo-Johnson Log-Likelihood Functions The scipy documention lists expressions for the Log-likelihood functions for the Box-Cox and Yeo-Johnson transformations here and here. I'm looking for a source that explains how one arrives at these expressions. I would like to experiment with my own transformations, and am therefore keen to understand how one goes from a data transformation to a Log-likelihood function.

Box-Cox Transformation: Parametric family of transformations $y\mapsto y^{(\lambda) }$ defined by
\begin{align}y^{(\lambda)} &:=\begin{cases}\frac{y^\lambda-1}{\lambda}~~&\lambda \ne 0\\ \ln y~~&\lambda =0\end{cases}\tag{1.I}\\y^{(\lambda)} &:=\begin{cases}\frac{(y+\lambda_2)^{\lambda_1}-1}{\lambda_1}~~&\lambda_1 \ne 0\\ \ln( y+\lambda_2)~~&\lambda_1=0;\end{cases}\tag{1.II}\end{align}
the former holds for $y>0$ and the latter for $y>-\lambda_2.$
Now $\mathbf y^{(\lambda) }\sim\mathcal N(\mathbf X\boldsymbol\beta, \sigma^2\mathbf I). $ Let $\boldsymbol\theta:=(\lambda, \beta, \sigma^2). $ Therefore, for fixed $\lambda,$ \begin{align}\mathcal L(\boldsymbol\theta|\mathbf y, \mathbf X) &=\frac{\exp{\left[-\frac{\left(\mathbf y^{(\lambda)}-\mathbf X\boldsymbol\beta\right)^\mathsf T\left(\mathbf y^{(\lambda)}-\mathbf X\boldsymbol\beta\right)}{2\sigma^2}\right]}}{(2\pi\sigma^2)^{n/2}}~\mathcal J(\lambda;\mathbf y)\tag{2.I}\\\implies \ln \mathcal L_\max&=C -\frac n2\ln \hat{\sigma}^2+ \ln \mathcal J,\tag {2.II}\end{align}
where $\mathcal J:= \prod \left|\frac{\mathrm d}{\mathrm dy_i}y_i^{(\lambda)}\right|:$ for $\mathrm{ (1.I)}, ~\ln \mathcal J= (\lambda-1)\sum\ln y_i$ and for $\mathrm{(1.II) }~\ln \mathcal J = (\lambda_1-1)\sum\ln( y_i+\lambda_2) ;~\hat{\sigma}^2= \frac{{\mathbf y^{(\lambda)}}^\mathsf T\left[\mathbf I-\mathbf X\left(\mathbf X^\mathsf T\mathbf X\right) ^{-}\mathbf X\right]\mathbf y^{(\lambda)}}{n}.$
Yeo-Johnson Transformation: $ y\mapsto \psi(\lambda, y) $ as
\begin{align}\psi(y,\lambda)&:= \begin{cases}\frac{(y+1)^\lambda-1}{\lambda}~~&(y\geq 0,~\lambda \ne 0) \\\ln(y+1)~&  (y\geq 0,~\lambda = 0) \\ -\frac{(-y+1)^{2-\lambda}-1}{2-\lambda}~~&(y< 0,~\lambda \ne 2)\\ \ln(-y+1) ~&  (y< 0,~\lambda = 2)\end{cases}.\tag 3\end{align}
By assumption, the transformed variables $\{\psi(\lambda,y_i)\}$ are normally distributed with mean, say, $\mu$ and variance $\sigma^2.$ Therefore, akin to $(2) ,$ $$\ln \mathcal L= -\frac{ n}{2  }\ln{2\pi} -\frac{n}{2}\ln \sigma^2-\frac{ 1}{2\sigma^2}\sum(\psi(\lambda, y_i) -\mu)^2+(\lambda-1) \sum\operatorname{sgn}(y_i) \ln(|y_i|+1). \tag 4$$ The corresponding profile loglikelihood would then be maximised using the MLE of $\mu, ~\sigma^2$ which are routine procedures away to compute in same way as above.

References:
$\rm [I]$ An Analysis of Transformations, G. E. P. Box, D. R. Cox, Journal of the Royal Statistical Society. Series B (Methodological) $26,$ no. $2 ~(1964): 211–52;$ pp. $214-216.$
$\rm[II]$ A New Family of Power Transformations to Improve Normality or Symmetry, In-Kwon Yeo, Richard A. Johnson, Biometrika $87,$ no. $4 ~(2000): 954–59;$ pp. $956-958.$
