Replies: 8 comments 22 replies
-
|
@tvwenger A PR would be appreciated! |
Beta Was this translation helpful? Give feedback.
-
|
|
Beta Was this translation helpful? Give feedback.
-
|
|
Beta Was this translation helpful? Give feedback.
-
|
@tvwenger from reading this is just a scaled chi distribution with 3 degrees of freedom? If that's the case you don't need to implement a new standard distribution and should instead use a import pymc as pm
def maxwell_dist(scale, size):
return pm.math.sqrt(pm.ChiSquared.dist(nu=3, size=size)) * scale
with pm.Model():
...
pm.CustomDist("x", scale, dist=maxwell_dist, **kwargs)If you want to implement it in PyMC-Experimental you can make a helper class that returns such CustomDist, something like class Maxwell:
def __new__(name, scale, **kwargs):
if "observed" not in kwargs:
kwargs.set_default("transform", pm.distributions.transforms.log)
return CustomDist(name, scale, dist=maxwell_dist, class_name="Maxwell", **kwargs)
def dist(scale, **kwargs):
return CustomDist.dist(scale, dist=maxwell_dist, class_name="Maxwell", **kwargs)Testing will also be easier, just check the logp matches the same as scipy in one or two points. |
Beta Was this translation helpful? Give feedback.
-
|
@tvwenger That's why there is a |
Beta Was this translation helpful? Give feedback.
-
|
@twiecki it's not just a reparam of an RV, it's a transformation of one, and the suggestion of using a |
Beta Was this translation helpful? Give feedback.
-
|
Here's a demonstration: import numpy as np
from scipy.stats import chi, chi2
import matplotlib.pyplot as plt
x = np.linspace(0.0, 5.0, 100)
pdf_chi = chi(3, loc=0.0, scale=1.0).pdf(x)
pdf_chi2 = chi2(3, loc=0.0, scale=1.0).pdf(x)
fig, ax = plt.subplots()
ax.plot(x, pdf_chi, 'k-', label=r'$\chi$')
ax.plot(x, pdf_chi2, 'r-', label=r'$\chi^2$')
ax.plot(x, np.sqrt(pdf_chi2), 'b-', label=r'$\sqrt{\chi^2}$')
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
ax.legend(loc='best')
fig.tight_layout()
fig.show()
plt.close(fig)I'm working on a PR to implement the chi distribution in |
Beta Was this translation helpful? Give feedback.
-
|
PyMC CustomDist is not taking the sqrt of the pdf of chi2. It is doing the right change of variables, taking the square root of the values and adding the jacobian of the transformation. Here is a visual proof that PyMC does the right thing: import scipy.stats as st
import pymc as pm
import numpy as np
test_vals = np.linspace(0, 5)
def maxwell_dist(scale, size):
return pm.math.sqrt(pm.ChiSquared.dist(nu=3, size=size)) * scale
scale = 1
plt.plot(pm.logp(pm.CustomDist.dist(scale, dist=maxwell_dist), test_vals).exp().eval(), label="PyMC CustomDist")
# multiply by 1.01 so we can see there are actually two lines
plt.plot(st.maxwell(scale=scale).pdf(test_vals) *1.01, label="Scipy Maxwell")
plt.legend()
|
Beta Was this translation helpful? Give feedback.
-
|
Great, thanks for the clarification! |
Beta Was this translation helpful? Give feedback.
-
|
@ricardoV94 I took a stab at an implementation following your suggestion in |
Beta Was this translation helpful? Give feedback.
-
|
Awesome, I'll take a look soon. |
Beta Was this translation helpful? Give feedback.
-
|
@ricardoV94 I'm starting to wonder if a |
Beta Was this translation helpful? Give feedback.
-
|
@ricardoV94 The problem continues! I thought the issue was that the import pymc as pm
print(pm.__version__)
import numpy as np
def func(obs, x):
return x * pm.math.sin(obs)
def main():
obs1 = np.deg2rad(np.array([20.0, 30.0, 50.0]))
obs2 = np.array([1.0, 2.0, 3.0])
with pm.Model() as model:
scale = pm.HalfNormal("scale", sigma=1.0)
def maxwell_dist(scale, size):
return pm.math.sqrt(pm.ChiSquared.dist(nu=3, size=size)) * scale
x = pm.CustomDist("x", scale, dist=maxwell_dist, shape=(3,), initval=np.ones(3))
mu = func(obs1, x)
pm.Normal("obs", mu=mu, sigma=1.0, observed=obs2)
with model:
trace = pm.sample()
main() |
Beta Was this translation helpful? Give feedback.
-
|
@tvwenger that's a great bug find! I am opening a PR to solve it |
Beta Was this translation helpful? Give feedback.
-
|
@ricardoV94 Thanks for submitting the PR! Any simple workaround until it's fixed in dev? |
Beta Was this translation helpful? Give feedback.
-
|
Yes, you can pass a custom logp function to CustomDist that looks like this: def maxwell_logp(value, *params):
# recreate RV
rv = maxwell_dist(*params, size=value.shape)
# Ask PyMC to derive its logp
return pm.logp(rv, value)
...
x = CustomDist(..., logp=maxwell_logp)I used the Maxwell example, but should work for whatever CustomDist you are playing with |
Beta Was this translation helpful? Give feedback.
-
|
@ricardoV94 Thanks for the suggestion! Since I really need a truncated Maxwell distribution for my model, I guess I also have to pass a import pymc as pm
print(pm.__version__)
import numpy as np
from scipy.stats import chi
def maxwell_dist(scale, size=None):
return pm.math.sqrt(pm.ChiSquared.dist(nu=3, size=size)) * scale
def maxwell_logp(value, *params, **kwargs):
rv = maxwell_dist(*params, size=value.shape)
return pm.logp(rv, value)
def maxwell_random(scale, rng=None, size=None):
return chi.rvs(3, loc=0, scale=scale, size=size, random_state=rng)
def main():
scale_true = 2.0
a = 5.0
rng = np.random.RandomState(seed=1234)
obs = chi.rvs(3, scale=scale_true/a, size=10, random_state=rng)
with pm.Model() as model:
scale = pm.HalfNormal("scale", sigma=1.0)
x = pm.Truncated(
"x",
pm.CustomDist.dist(
scale/a,
logp=maxwell_logp,
random=maxwell_random,
),
upper=5.0,
observed=obs,
)
with model:
trace = pm.sample(init="adapt_diag")
main()Fails with I've attempted to implement the import pymc as pm
print(pm.__version__)
import numpy as np
from scipy.stats import chi
def maxwell_dist(scale, size=None):
return pm.math.sqrt(pm.ChiSquared.dist(nu=3, size=size)) * scale
def maxwell_logp(value, *params, **kwargs):
rv = maxwell_dist(*params, size=value.shape)
return pm.logp(rv, value)
def maxwell_random(scale, rng=None, size=None):
return chi.rvs(3, loc=0, scale=scale, size=size, random_state=rng)
def maxwell_logcdf(value, *params, **kwargs):
rv = maxwell_dist(*params, size=value.shape)
return pm.logcdf(rv, value)
def main():
scale_true = 2.0
# a = 5.0 # this works always
a = 1.0 + np.arange(10) # this fails with Truncated CustomDist
rng = np.random.RandomState(seed=1234)
obs = chi.rvs(3, scale=scale_true/a, size=10, random_state=rng)
with pm.Model() as model:
scale = pm.HalfNormal("scale", sigma=1.0)
x = pm.Truncated(
"x",
pm.CustomDist.dist(
scale/a,
logp=maxwell_logp,
random=maxwell_random,
logcdf=maxwell_logcdf,
),
upper=5.0,
observed=obs,
)
with model:
trace = pm.sample(init="adapt_diag")
main()This fails with an obscure error if |
Beta Was this translation helpful? Give feedback.
-
|
Implemented in pymc-devs/pymc-extras#261 |
Beta Was this translation helpful? Give feedback.
-
I have need of a Maxwell distribution in
pymc. I don't think any of the existing distributions work for modelling a Maxwell distribution. Unless there is already an implementation, I can work on a PR for it.Beta Was this translation helpful? Give feedback.
All reactions