Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Symbolic dimensions allowed in MultivariateEwens #16914

Merged
merged 3 commits into from Jun 1, 2019
Merged
Changes from 2 commits
Commits
File filter...
Filter file types
Jump to…
Jump to file or symbol
Failed to load files and symbols.

Always

Just for now

@@ -16,6 +16,7 @@
from sympy import (Basic, Lambda, sympify, Indexed, Symbol, ProductSet, S,
Dummy)
from sympy.concrete.summations import Sum, summation
from sympy.concrete.products import Product
from sympy.core.compatibility import string_types
from sympy.core.containers import Tuple
from sympy.integrals.integrals import Integral, integrate
@@ -64,7 +65,11 @@ def value(self):
@property
def component_count(self):
_set = self.distribution.set
return len(_set.args) if isinstance(_set, ProductSet) else 1
if isinstance(_set, ProductSet):
return S(len(_set.args))
elif isinstance(_set, Product):
return _set.limits[0][-1]
return S(1)

@property
def pdf(self):
@@ -87,6 +92,9 @@ def symbols(self):

def marginal_distribution(self, *indices):
count = self.component_count
if count.atoms(Symbol):
raise ValueError("Marginal distributions cannot be computed "
"for symbolic dimensions. It is a work under progress.")
orig = [Indexed(self.symbol, i) for i in range(count)]
all_syms = [Symbol(str(i)) for i in orig]
replace_dict = dict(zip(all_syms, orig))
@@ -182,7 +190,7 @@ class JointRandomSymbol(RandomSymbol):
"""
def __getitem__(self, key):
if isinstance(self.pspace, JointPSpace):
if self.pspace.component_count <= key:
if (self.pspace.component_count <= key) == True:
raise ValueError("Index keys for %s can only up to %s." %
(self.name, self.pspace.component_count - 1))
return Indexed(self, key)
@@ -1,7 +1,7 @@
from sympy import (sympify, S, pi, sqrt, exp, Lambda, Indexed, Gt, IndexedBase,
besselk, gamma, Interval, Range, factorial, Mul, Integer,
Add, rf, Eq, Piecewise, ones, Symbol, Pow, Rational, Sum,
imageset, Intersection, Matrix)
imageset, Intersection, Matrix, symbols, Product, IndexedBase)
from sympy.matrices import ImmutableMatrix
from sympy.matrices.expressions.determinant import det
from sympy.stats.joint_rv import (JointDistribution, JointPSpace,
@@ -339,25 +339,41 @@ class MultivariateEwensDistribution(JointDistribution):

@staticmethod
def check(n, theta):
_value_check(isinstance(n, Integer) and (n > 0) == True,
_value_check((n > 0),
"sample size should be positive integer.")
_value_check(theta.is_positive, "mutation rate should be positive.")

@property
def set(self):
prod_set = Range(0, self.n//1 + 1)
if not isinstance(self.n, Integer):
i = Symbol('i', integer=True, positive=True)
return Product(Intersection(S.Naturals0, Interval(0, self.n//i)),
(i, 1, self.n))
prod_set = Range(0, self.n + 1)
for i in range(2, self.n + 1):
prod_set *= Range(0, self.n//i + 1)
return prod_set

def pdf(self, *syms):
n, theta = self.n, self.theta
condi = isinstance(self.n, Integer)
if not (isinstance(syms[0], IndexedBase) or condi):
raise ValueError("Please use IndexedBase object for syms as "
"the dimension is symbolic")
term_1 = factorial(n)/rf(theta, n)
term_2 = Mul.fromiter([theta**syms[j]/((j+1)**syms[j]*factorial(syms[j]))
for j in range(n)])
cond = Eq(sum([(k+1)*syms[k] for k in range(n)]), n)
if condi:
term_2 = Mul.fromiter([theta**syms[j]/((j+1)**syms[j]*factorial(syms[j]))
for j in range(n)])
cond = Eq(sum([(k + 1)*syms[k] for k in range(n)]), n)
return Piecewise((term_1 * term_2, cond), (0, True))
syms = syms[0]
j, k = symbols('j, k', positive=True, integer=True)
term_2 = Product(theta**syms[j]/((j+1)**syms[j]*factorial(syms[j])),
(j, 0, n - 1))
cond = Eq(Sum((k + 1)*syms[k], (k, 0, n - 1)), n)
return Piecewise((term_1 * term_2, cond), (0, True))


def MultivariateEwens(syms, n, theta):
"""
Creates a discrete random variable with Multivariate Ewens
@@ -1,5 +1,6 @@
from sympy import (symbols, pi, oo, S, exp, sqrt, besselk, Indexed, Sum, simplify,
Mul, Rational, Integral, factorial, gamma, Piecewise, Eq)
Mul, Rational, Integral, factorial, gamma, Piecewise, Eq, Product,
IndexedBase)
from sympy.core.numbers import comp
from sympy.stats import density
from sympy.stats.joint_rv import marginal_distribution
@@ -66,7 +67,6 @@ def test_GeneralizedMultivariateLogGammaDistribution():
[h, h, h, 1]])
v, l, mu = (4, [1, 2, 3, 4], [1, 2, 3, 4])
y_1, y_2, y_3, y_4 = symbols('y_1:5', real=True)
n = symbols('n', negative=False, integer=True)
delta = symbols('d', positive=True)
G = GMVLGO('G', omega, v, l, mu)
Gd = GMVLG('Gd', delta, v, l, mu)
@@ -135,7 +135,10 @@ def test_MultivariateBeta():

def test_MultivariateEwens():
from sympy.stats.joint_rv_types import MultivariateEwens
n, theta = symbols('n theta', positive=True)

n, theta, i = symbols('n theta i', positive=True)

# tests for integer dimensions
theta_f = symbols('t_f', negative=True)
a = symbols('a_1:4', positive = True, integer = True)
ed = MultivariateEwens('E', 3, theta)
@@ -150,7 +153,14 @@ def test_MultivariateEwens():
(theta + 2)*factorial(a[1])),
Eq(2*a[1] + 1, 3)), (0, True))
raises(ValueError, lambda: MultivariateEwens('e1', 5, theta_f))
raises(ValueError, lambda: MultivariateEwens('e1', n, theta))

# tests for symbolic dimensions
eds = MultivariateEwens('E', n, theta)
a = IndexedBase('a')
den = ("Piecewise((factorial(n)*Product(theta**a[j]*(j + 1)**(-a[j])/"
"factorial(a[j]), (j, 0, n - 1))/RisingFactorial(theta, n),"
" Eq(n, Sum((k + 1)*a[k], (k, 0, n - 1)))), (0, True))")
assert str(density(eds)(a)) == den
This conversation was marked as resolved by czgdp1807

This comment has been minimized.

Copy link
@Upabjojr

Upabjojr May 31, 2019

Contributor

have you given a try to dummy_eq?

This comment has been minimized.

Copy link
@czgdp1807

czgdp1807 May 31, 2019

Author Member

Thanks it worked fine.


def test_Multinomial():
from sympy.stats.joint_rv_types import Multinomial
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.