custom stats distribution

[spl...@gmail.com]

Hello,

I'm trying to set up a custom distribution with the following function describing its PDF:

def pdf(x, beta, alpha, eta):
    t0 = beta / (alpha * (1 - x - eta))
    t1 = pow(np.log(1 - x - eta) / alpha, beta - 1)
    t2 = np.exp(-pow(np.log(1 - x - eta) / alpha, beta))
    return t0 * t1 * t2

It's a particular parameterization of the Log-Weibull distribution, with beta and alpha being the shape and scale parameters, respectively, and eta acts as a shift parameter representing the upper bound of the sample space, which is then defined as -inf < x < eta.  I'm interested in drawing random variates, computing the above PDF, and CDF.

I have limited experience with sympy, and none with sympy.stats, but it seems as if what I need is sympy.stats.ContinuousDistributionHandmade.  However, I'd appreciate any pointers.

Thanks,

--

Seb

[Gagandeep Singh (B17CS021)]

Hi, 

I think for creating RVs with custom distributions, one can use, 
ContinuousRV, DiscreteRV Or FiniteRV depending on the type of distribution. You can take a look at https://docs.sympy.org/latest/modules/stats.html#examples 

Note that the PDF/PMF should be a SymPy object. 

With Regards,
Gagandeep Singh
Github - https://github.com/czgdp1807
LinkedIn - https://www.linkedin.com/in/czgdp1807

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/8f73147b-7032-464d-b95c-a1b68537c1a7n%40googlegroups.com.

[spl...@gmail.com]

Here is my attempt with ContinuousRV:

from sympy.stats import ContinuousRV
from sympy import Lambda, log, Pow, exp, symbols, Interval, oo

# Parameters
_BETA = 10.56
_ETA = -90
_ALPHA = 3.11

x, beta, alpha, eta = symbols("x, beta, alpha, eta")

log_weibull_pdf = ((beta / (alpha * (1 - x - eta))) *
                   (Pow(log(1 - x - eta) / alpha, beta - 1)) *
                   (exp(-Pow(log(1 - x - eta) / alpha, beta))))
X = ContinuousRV(x, log_weibull_pdf.subs({beta: _BETA,
                                          alpha: _ALPHA,
                                          eta: _ETA}),
                 set=Interval(-oo, -_ETA))

which fails with:

ValueError: x**w where w is irrational is not defined for negative x

--

Seb

[Gagandeep Singh (B17CS021)]

IPython console for SymPy 1.7.dev (Python 3.6.9-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at https://docs.sympy.org/dev


In [1]: from sympy.stats import ContinuousRV

In [2]: from sympy import Lambda, log, Pow, exp, symbols, Interval, oo

In [3]: _BETA = 10.56

In [5]: _ETA = -90

In [6]: _ALPHA = 3.11

In [7]: x, beta, alpha, eta = symbols("x, beta, alpha, eta")

In [8]: log_weibull_pdf = ((beta / (alpha * (1 - x - eta))) *
   ...:                    (Pow(log(1 - x - eta) / alpha, beta - 1)) *
   ...:                    (exp(-Pow(log(1 - x - eta) / alpha, beta))))

In [9]: log_weibull_pdf.subs({beta: _BETA,
   ...:                                           alpha: _ALPHA,
   ...:                                           eta: _ETA}
   ...:
   ...:
   ...:
   ...:
   ...:                                           )
Out[9]:
                                             10.56                        
                     -6.25821230086363e-6⋅log     (91 - x)    9.56        
6.60867218971199e-5⋅ℯ                                     ⋅log    (91 - x)
──────────────────────────────────────────────────────────────────────────
                                  91 - x                                  

In [10]: pdf = _

In [11]: X = ContinuousRV(x, pdf, set=Interval(-oo, -_ETA))

In [12]: X
Out[12]: x

In [14]: from sympy.stats import density

In [15]: density(X)(y)
Out[15]:
⎧                                             10.56                          
⎪                     -6.25821230086363e-6⋅log     (91 - y)    9.56          
⎪6.60867218971199e-5⋅ℯ                                     ⋅log    (91 - y)  
⎨──────────────────────────────────────────────────────────────────────────  f
⎪                                  91 - y                                    
⎪                                                                            
⎩                                    0                                        

                 
                 
                 
or y ≥ -∞ ∧ y ≤ 90
                 
                 
    otherwise

Is the above what you were trying to do?     

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/16120812-5243-4f41-869e-7d03288dfb97n%40googlegroups.com.

--

With regards,

Gagandeep Singh

Github - https://github.com/czgdp1807/

Linkedin - https://www.linkedin.com/in/czgdp1807/

[Gagandeep Singh (B17CS021)]

Hi,

Try using the following,

X = ContinuousRV(x, pdf, set=Interval(-oo, -_ETA), check=False)

The keyword argument `check=True` in your case and hence the PDF is integrated for validation which fails. You can avoid that check by setting it to `False`. In the next release, it would be by default, `False`. See, https://github.com/sympy/sympy/pull/20120 and hence your code worked on the development version.  

[spl...@gmail.com]

It seems you're working with the dev version.  I'm using the latest release (1.6.2), where `ContinuousRV` does not have a `check` argument.  Why would the integration fail in the current release?

Thanks for the feedback,

--

Seb

[Gagandeep Singh (B17CS021)]

> Why would the integration fail in the current release?

It fails in master(dev version) too. It might be a problem with the integral module.

> I'm using the latest release (1.6.2), where `ContinuousRV` does not have a `check` argument.

I see, maybe `check` argument will also come with the next release. Since, we ensure CI checks pass before merge, you may try using `master`(dev version) directly. The following steps can be helpful for installing the dev version. I have tried them on Linux. See https://docs.sympy.org/dev/install.html for more details.

1. git clone https://github.com/sympy/sympy

2. cd /path/to/sympy/sympy

3. python3 setup.py develop

This way you will be able to do your task without much delay.

Thank you and apologies for inconveniences.  

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/4f6cd660-20e1-4fce-9f58-30ff3a65e505n%40googlegroups.com.