<p dir="auto">When trying to compute variance of <code class="notranslate">stats.truncnorm</code> it gives incorrect results when shape parameters <code class="notranslate">a</code> and <code class="notranslate">b</code> are arrays.</p>
<p dir="auto">It works correctly when <code class="notranslate">a</code> and <code class="notranslate">b</code> are floats.</p>
<p dir="auto">I am working with <code class="notranslate">scipy-1.6.1</code> and Python 3.9.1.</p>
<h4 dir="auto">Reproducing code example:</h4>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="import numpy as np
from scipy import stats

a = np.arange(10, dtype=float)
b = a + 1

# This will print correct mean
print('mean:', stats.truncnorm.stats(a, b, moments='m'))

# This will print incorrect variance
print('var:', stats.truncnorm.stats(a, b, moments='v'))

# This will print correct variance
print('var:', [stats.truncnorm.stats(aa, bb, moments='v') for (aa, bb) in zip(a, b)])"><pre class="notranslate"><span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>
<span class="pl-k">from</span> <span class="pl-s1">scipy</span> <span class="pl-k">import</span> <span class="pl-s1">stats</span>

<span class="pl-s1">a</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">arange</span>(<span class="pl-c1">10</span>, <span class="pl-s1">dtype</span><span class="pl-c1">=</span><span class="pl-s1">float</span>)
<span class="pl-s1">b</span> <span class="pl-c1">=</span> <span class="pl-s1">a</span> <span class="pl-c1">+</span> <span class="pl-c1">1</span>

<span class="pl-c"># This will print correct mean</span>
<span class="pl-en">print</span>(<span class="pl-s">'mean:'</span>, <span class="pl-s1">stats</span>.<span class="pl-s1">truncnorm</span>.<span class="pl-en">stats</span>(<span class="pl-s1">a</span>, <span class="pl-s1">b</span>, <span class="pl-s1">moments</span><span class="pl-c1">=</span><span class="pl-s">'m'</span>))

<span class="pl-c"># This will print incorrect variance</span>
<span class="pl-en">print</span>(<span class="pl-s">'var:'</span>, <span class="pl-s1">stats</span>.<span class="pl-s1">truncnorm</span>.<span class="pl-en">stats</span>(<span class="pl-s1">a</span>, <span class="pl-s1">b</span>, <span class="pl-s1">moments</span><span class="pl-c1">=</span><span class="pl-s">'v'</span>))

<span class="pl-c"># This will print correct variance</span>
<span class="pl-en">print</span>(<span class="pl-s">'var:'</span>, [<span class="pl-s1">stats</span>.<span class="pl-s1">truncnorm</span>.<span class="pl-en">stats</span>(<span class="pl-s1">aa</span>, <span class="pl-s1">bb</span>, <span class="pl-s1">moments</span><span class="pl-c1">=</span><span class="pl-s">'v'</span>) <span class="pl-k">for</span> (<span class="pl-s1">aa</span>, <span class="pl-s1">bb</span>) <span class="pl-c1">in</span> <span class="pl-en">zip</span>(<span class="pl-s1">a</span>, <span class="pl-s1">b</span>)])</pre></div>
<p dir="auto"><strong>Output of this program:</strong></p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="mean: array([0.45986223, 1.38316905, 2.31582133, 3.26045429, 4.21683078,
       5.18314709, 6.1572109 , 7.13706716, 8.12118899, 9.10845629])
var: array([-8.60292538, -8.60292538, -8.60292538, -8.60292538, -8.60292538,
       -8.60292538, -8.60292538, -8.60292538, -8.60292538, -8.60292538])
var: [array(0.07965182), array(0.07274289), array(0.06152078), array(0.04927779), array(0.03829029), array(0.02945243), array(0.02274838), array(0.01779288), array(0.01414854), array(0.01144881)]"><pre class="notranslate"><code class="notranslate">mean: array([0.45986223, 1.38316905, 2.31582133, 3.26045429, 4.21683078,
       5.18314709, 6.1572109 , 7.13706716, 8.12118899, 9.10845629])
var: array([-8.60292538, -8.60292538, -8.60292538, -8.60292538, -8.60292538,
       -8.60292538, -8.60292538, -8.60292538, -8.60292538, -8.60292538])
var: [array(0.07965182), array(0.07274289), array(0.06152078), array(0.04927779), array(0.03829029), array(0.02945243), array(0.02274838), array(0.01779288), array(0.01414854), array(0.01144881)]
</code></pre></div>
<p dir="auto"><code class="notranslate">truncnorm.mean</code> and <code class="notranslate">truncnorm.var</code> appear to give incorrect outputs when <code class="notranslate">loc</code> or <code class="notranslate">scale</code> are array-like and contain more than one number:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; import scipy.stats as stats
&gt;&gt;&gt; stats.truncnorm.mean(-1, 2, loc=0)
0.229637179091329
&gt;&gt;&gt; stats.truncnorm.mean(-1, 2, loc=[0])
array([0.22963718])
&gt;&gt;&gt; stats.truncnorm.mean(-1, 2, loc=[0, 0]) # expect to get 0.2296 twice
array([0., 0.])
&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=0)
0.29112509477279314
&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=[0])
array([0.29112509])
&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=[0, 0]) # expect to get 0.2911 twice
array([1., 1.])"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; import scipy.stats as stats
&gt;&gt;&gt; stats.truncnorm.mean(-1, 2, loc=0)
0.229637179091329
&gt;&gt;&gt; stats.truncnorm.mean(-1, 2, loc=[0])
array([0.22963718])
&gt;&gt;&gt; stats.truncnorm.mean(-1, 2, loc=[0, 0]) # expect to get 0.2296 twice
array([0., 0.])
&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=0)
0.29112509477279314
&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=[0])
array([0.29112509])
&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=[0, 0]) # expect to get 0.2911 twice
array([1., 1.])
</code></pre></div>
<p dir="auto"><code class="notranslate">truncnorm.moment</code> doesn't seem to support array-likes at all, either:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; stats.truncnorm.moment(1, -1, 1, loc=[0])
Traceback (most recent call last):
  File &quot;&lt;stdin&gt;&quot;, line 1, in &lt;module&gt;
  File &quot;C:\ProgramData\Miniconda3\lib\site-packages\scipy\stats\_distn_infrastructure.py&quot;, line 1266, in moment
    fac = float(scale) / float(loc)
TypeError: float() argument must be a string or a number, not 'list'"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; stats.truncnorm.moment(1, -1, 1, loc=[0])
Traceback (most recent call last):
  File "&lt;stdin&gt;", line 1, in &lt;module&gt;
  File "C:\ProgramData\Miniconda3\lib\site-packages\scipy\stats\_distn_infrastructure.py", line 1266, in moment
    fac = float(scale) / float(loc)
TypeError: float() argument must be a string or a number, not 'list'
</code></pre></div>
<p dir="auto">This is on Scipy version 1.6.0. This problem was even worse on version 1.4.1, where <code class="notranslate">truncnorm.var()</code> could produce negative numbers:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=0)
0.29112509477279314
&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=[0, 0])
array([-0.41774981, -0.41774981])"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=0)
0.29112509477279314
&gt;&gt;&gt; stats.truncnorm.var(-1, 1, loc=[0, 0])
array([-0.41774981, -0.41774981])
</code></pre></div>
<p dir="auto">I don't think I'm using these functions incorrectly, as the docstrings for <code class="notranslate">mean()</code>, <code class="notranslate">var()</code>, and <code class="notranslate">moment()</code> all say that <code class="notranslate">loc</code> and <code class="notranslate">scale</code> can be array-like.</p>