Using scipy.stats library or another method to generate data follows a distribution in a specific boundary
Question: I want to sample with `scipy.stats` library, using an upper and a lower
boundary for the sampled data. I am interested to use `scipy.stats.lognorm`
and `scipy.stats.expon` and set a constrain `(low<=x<=up)` on the limits of
generated data points and also estimate `logp` with considering these limits.
For instance, I can not do

    
    
    LogNormal=scipy.stats.lognorm(q=[0,5],scale=[0.25],loc=0.0) #q:upper and lower limits, scale=sigma, loc=mean
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "/vol/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py", line 739, in __call__
        return self.freeze(*args, **kwds)
      File "/vol/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py", line 736, in freeze
        return rv_frozen(self, *args, **kwds)
      File "/vol/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py", line 434, in __init__
        shapes, _, _ = self.dist._parse_args(*args, **kwds)
    TypeError: _parse_args() got an unexpected keyword argument 'q'
    

The documentation is a bit confusing, which one is `sigma` and which input
parameter is `mean`? Could anybody give an example, how they should be set
with boundaries?


Answer: There are several problems in your implementation

1, your pdf can not be evaluated at x=0

2,
`-log(1./sqrt(2*pi)/self.sigma*exp(-0.5*((log(value)-self.mu)/self.sigma)**2))`
should be:
`-log(1./sqrt(2*pi)/self.sigma/value*exp(-0.5*((log(value)-self.mu)/self.sigma)**2))`

(And there may be more)

Another consideration is that you may want to keep the parameterization the
same as `scipy` to avoid future confusion.

Therefore, a minimal implementation:

    
    
    In [112]:
    import scipy.stats as ss
    import scipy.optimize as so
    import numpy as np
    
    class bounded_distr(object):
        def __init__(self, parent_dist):
            self.parent = parent_dist
        def bnd_lpdf(self, x, limits=None, *args, **kwargs):
            if limits and np.diff(limits)<=0:
                return -np.inf #nan may be better idea
            else:
                _v = -log(self.parent.pdf(x, *args, **kwargs))
                _v[x<=limits[0]] = -np.inf
                _v[x>=limits[1]] = -np.inf
                return _v
        def bnd_cdf(self, x, limits=None, *args, **kwargs):
            if limits and np.diff(limits)<=0:
                return 0 #nan may be better idea
            elif limits:
                _v1 = self.parent.cdf(x, *args, **kwargs)
                _v2 = self.parent.cdf(limits[0], *args, **kwargs)
                _v3 = self.parent.cdf(limits[1], *args, **kwargs)
                _v4 = (_v1-_v2)/(_v3-_v2)
                _v4[_v4<0] = np.nan
                _v4[_v4>1] = np.nan
                return _v4
            else:
                return self.parent.cdf(x, *args, **kwargs)
        def bnd_rvs(self, size, limits=None, *args, **kwargs):
            if limits and np.diff(limits)<=0:
                return np.repeat(np.nan, size) #nan may be better idea
            elif limits:
                low, high = limits
                rnd_cdf = np.random.uniform(self.parent.cdf(x=low, *args, **kwargs),
                                            self.parent.cdf(x=high, *args, **kwargs),
                                            size=size)
                return self.parent.ppf(q=rnd_cdf, *args, **kwargs)
            else:
                return self.parent.rvs(size=size, *args, **kwargs)
    In [113]:
    
    bnd_logn = bounded_distr(ss.lognorm)
    In [114]:
    
    bnd_logn.bnd_rvs(10, limits=(0.1, 0.9), s=1, loc=0)
    Out[114]:
    array([ 0.23167598,  0.43185726,  0.34763109,  0.71020467,  0.5216074 ,
            0.60883528,  0.34353607,  0.84530444,  0.64145739,  0.82082447])
    In [115]:
    
    bnd_logn.bnd_lpdf(np.linspace(0,1,10), limits=(0.1, 0.9), s=1, loc=0)
    Out[115]:
    array([        inf,  1.13561188,  0.54598554,  0.42380072,  0.43681222,
            0.50389845,  0.5956744 ,  0.69920358,  0.80809192,  0.91893853])
    In [116]:
    
    bnd_logn.bnd_cdf(np.linspace(0,1,10), limits=(0.1, 0.9), s=1, loc=0)
    Out[116]:
    array([        nan,  0.00749028,  0.12434152,  0.28010562,  0.44267888,
            0.59832448,  0.74188947,  0.87201574,  0.98899161,         nan])
    

