Transforming a general distribution to that of a power

You want to call CustomDist.dist from the dist method, and not use the dist kwarg…

However you should actually use it. You get everything for free by implementing a CustomDist with a dist kwarg instead of with a logp and random kwargs. Probably sounds confusing so let me show with one less layer of abstraction.

(untested code)

import pymc as pm

def square_invgamma_dist(mu, sigma, size):
  return pm.InverseGamma.dist(mu=mu, sigma=sigma, size=size) ** 2

with pm.Model() as m:

  mu = 2
  sigma = 1
  sd_dist = pm.CustomDist.dist(mu, sigma, dist=square_invgamma_dist)
  pm.LKJCholeskyCov("lkj", n=3, eta=1, sd_dist=sd_didt)

The random method is trivially replaced because dist is already a symbolic random object. PyMC can then use this to infer the logp with the change of variables for the squared inversegamma automatically.

That relies on pymc logp inference which you can test directly by calling:

pm.logp(
  pm.InverseGamma.dist(mu=1, sigma=2) ** 2,
  3.5,
).eval()