<p>The easiest way is probably to leverage the truncated normal distribution as provided by Scipy.</p>
<p>This gives the following code, with ν (nu) as the variable of the standard Gaussian distribution, and τ (tau) mapping to ν<sub>0</sub> on that distribution. This function returns a Numpy array containing ranCount lognormal variates:</p>
<pre><code>import  numpy  as  np
from scipy.stats import truncnorm

def getMySamplesScipy(ranCount, mu, sigma, tau):
    nu0 = (math.log(tau) - mu) / sigma              # position of tau on unit Gaussian
    xs = truncnorm.rvs(nu0, np.inf, size=ranCount)  # truncated unit normal samples
    ys = np.exp(mu + sigma * xs)                    # go back to x space
    return ys
</code></pre>
<p>If for some reason this is not suitable, well some of the tricks commonly used for Gaussian variates, such as <a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform" rel="nofollow noreferrer">Box-Muller</a> do not work for a truncated distribution, but we can resort always to a general principle: the <a href="https://en.wikipedia.org/wiki/Inverse_transform_sampling" rel="nofollow noreferrer">Inverse Transform Sampling</a> theorem.</p>
<p>So we generate cumulative probabilities for our variates, by transforming uniform variates. And we trust Scipy, using its inverse of the <em>erf</em> error function to go back from our probabilities to the x space values.</p>
<p>This gives something like the following Python code (without any attempt at optimization):</p>
<pre><code>import  math
import  random
import  numpy  as  np
import  numpy.random  as  nprd
import  scipy.special as  spfn

# using the &quot;Inverse Method&quot;:
def getMySamples(ranCount, mu, sigma, tau):
    nu0 = (math.log(tau) - mu) / sigma      # position of tau in standard Gaussian curve
    headCP = (1/2) * (1 + spfn.erf(nu0/math.sqrt(2)))
    tailCP = 1.0 - headCP                   # probability of being in the &quot;tail&quot;

    uvs = np.random.uniform(0.0, 1.0, ranCount)  # uniform variates
    cps = (headCP + uvs * tailCP)                # Cumulative ProbabilitieS
    nus = (math.sqrt(2)) * spfn.erfinv(2*cps-1)  # positions in standard Gaussian
    xs  = np.exp(mu + sigma * nus)               # go back to x space
    return xs
</code></pre>
<h2 id="alternatives">Alternatives:</h2>
<p>We can leverage the significant amount of material related to the <a href="https://en.wikipedia.org/wiki/Truncated_normal_distribution" rel="nofollow noreferrer">Truncated Gaussian distribution</a>.</p>
<p>There is a relatively recent (2016) <a href="https://www.iro.umontreal.ca/%7Elecuyer/myftp/papers/vt16truncnormal.pdf" rel="nofollow noreferrer">review paper</a> on the subject by Zdravko Botev and Pierre L'Ecuyer. This paper provides a pointer to publicly available <a href="https://cran.r-project.org/web/packages/TruncatedNormal" rel="nofollow noreferrer">R source code</a>. Some material is seriously old, for example the 1986 book by Luc Devroye: <a href="http://www.eirene.de/Devroye.pdf" rel="nofollow noreferrer">Non-Uniform Random Variate Generation</a>.</p>
<p>For example, a possible rejection-based method: if τ (tau) maps to ν<sub>0</sub> on the standard Gaussian curve, the unit Gaussian distribution is like exp(-ν<sup>2</sup>/2). If we write ν = ν<sub>0</sub> + δ, this is proportional to: exp(-δ<sup>2</sup>/2) * exp(-ν<sub>0</sub>*δ).</p>
<p>The idea is to approximate the exact distribution beyond ν<sub>0</sub> by an <strong>exponential one</strong>, of parameter ν<sub>0</sub>. Note that the exact distribution is constantly <em>below</em> the approximate one. Then we can randomly accept the relatively cheap exponential variates with a probability of exp(-δ<sup>2</sup>/2).</p>
<p>We can just pick an equivalent algorithm in the literature. In the Devroye book, chapter IX page 382, there is some pseudo-code:</p>
<p>REPEAT
generate independent exponential random variates X and Y
UNTIL  X<sup>2</sup> &lt;= 2*ν<sub>0</sub><sup>2</sup>*Y</p>
<p>RETURN  R &lt;-- ν<sub>0</sub> + X/ν<sub>0</sub></p>
<p>for which a Numpy rendition could be written like this:</p>
<pre><code>def getMySamplesXpRj(rawRanCount, mu, sigma, tau):
    nu0 = (math.log(tau) - mu) / sigma      # position of tau in standard Gaussian
    if (nu0 &lt;= 0):
        print(&quot;Error: τ (tau) too small in getMySamplesXpRj&quot;)
    rnu0 = 1.0 / nu0

    xs  = nprd.exponential(1.0,  rawRanCount)   # exponential &quot;raw&quot; variates
    ys  = nprd.exponential(1.0,  rawRanCount)

    allSamples = nu0 + (rnu0 * xs)
    boolArray  = (xs*xs - 2*nu0*nu0*ys) &lt;= 0.0
    samples    = allSamples[boolArray]

    ys  = np.exp(mu + sigma * samples)               # go back to x space
    return ys
</code></pre>
<p>According to Table 3 in the Botev-L'Ecuyer paper, the rejection rate of this algorithm is nicely low.</p>
<p>Besides, if you are willing to allow for some sophistication, there is also some literature about the <a href="https://en.wikipedia.org/wiki/Ziggurat_algorithm" rel="nofollow noreferrer"><em>Ziggurat</em> algorithm</a> as used for truncated Gaussian distributions, for example the 2012 <a href="https://arxiv.org/pdf/1201.6140" rel="nofollow noreferrer">arXiv 1201.6140 paper</a> by Nicolas Chopin at ENSAE-CREST.</p>
<p><strong>Side note:</strong> with recent versions of Python, it seems that you can use Greek letters for your variable names directly, σ instead of sigma, τ instead of tau, just as in the statistics books:</p>
<pre><code>$ python3
Python 3.9.6 (default, Jun 29 2021, 00:00:00)
&gt;&gt;&gt; 
&gt;&gt;&gt; σ = 2
&gt;&gt;&gt; τ = 7
&gt;&gt;&gt; 
&gt;&gt;&gt; στ = σ * τ
&gt;&gt;&gt; 
&gt;&gt;&gt; στ + 1
15
&gt;&gt;&gt; 

</code></pre>
Scipy: sample lognormal distribution from tail only
<p>I'm building a simulation which requires random draws from the tail of a lognormal distribution. A threshold τ (tau) is chosen, and a resulting conditional distribution is given by:
<a href="https://i.stack.imgur.com/54uNl.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/54uNl.png" alt="Formula for truncated lognormal distribution" /></a></p>
<p>I need to randomly sample from that conditional distribution, where F(x) is lognormal with a chosen µ (mu) and σ (sigma), and τ (tau) is set by the user.</p>
<p>My inelegant solution right now is simply to sample from the lognormal, tossing out any values under τ (tau), until I have the sample size I need. But I'm sure this can be improved.</p>
<p>Thanks for the help!</p>