1

The PDF I have is:

$$ P(\varepsilon) = \frac{1}{4k} \mathrm{sech}^2 \frac{\varepsilon}{2k}$$

CDF:

$$ CDF(\varepsilon) = \frac{\tanh(\frac{\varepsilon}{2k})}{2} $$

and Inverse CDF:

$$ CDF^{-1}(c) = \frac{c}{2\tanh^{-1}(2k_BT)} $$

Now what should I do to generate N samples, with my PDF P(\varepsilon) of given k?

Here is what I have so far: 

import numpy as np
def sampleConduction(N, kT):

    for i in range(N):
        c = np.random.random()
        list = [c]
        print('list=',list) 
    g = np.column_stack(list)
    print('g=',g)

     d = np.arctanh(2*list)+2*kT

    return d

sampleConduction(3,0.1)

The result of this code is:

list= [0.7687170402304889]
list= [0.24759083266582882]
list= [0.18334770166799108]
g= [[0.1833477]]

array([0.38544466, 0.38544466]) # d

I wanted to combine the 3 lists into one, but failed by using np.column_stack(). What should I try next?

Thanks a lot!

g_tenga
  • 51
  • 7
  • Two comments, not related to the actual question that you ask: (1) Your formula for the inverse CDF does not look correct. That should be `*`, not `+`, in the formula. (2) The expression `tanh^{-1}(x)` is not `1/tanh(x)`. `tanh^{-1}` is *the function that is the inverse of the hyperbolic tangent function*. In numpy, this function is [`numpy.arctanh`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.arctanh.html). – Warren Weckesser Sep 16 '18 at 00:16
  • @WarrenWeckesser Thank you! I have corrected my math error. – g_tenga Sep 16 '18 at 00:40
  • Subclassing `scipy.stats.rv_continuous` may be an option. See e.g. [here](https://stackoverflow.com/a/22457754/7207392) for a starting point. This class can generate samples when given only the pdf, but I suspect it will perform better if you also provide the cdf or better still the inverse cdf. – Paul Panzer Sep 16 '18 at 02:11
  • FYI: Your distribution is known as the [logistic distribution](https://en.wikipedia.org/wiki/Logistic_distribution). It is implemented in scipy as [scipy.stats.logistic](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.logistic.html). – Warren Weckesser Sep 16 '18 at 04:24
  • If all you need are random variates, numpy has the function [`numpy.random.logistic`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.logistic.html). – Warren Weckesser Sep 16 '18 at 04:27
  • @WarrenWeckesser Thank you! numpy.random.logistic is indeed what I was looking for. I guess I have complicated this problem. – g_tenga Sep 16 '18 at 05:56

0 Answers0