Utilisation de la diffusion numérique avec scipy truncnorm

Question

Je souhaite évaluer une distribution normale tronquée unilatérale pour différentes valeurs du quantile et différentes valeurs de la moyenne non tronquée. Pour plus d'efficacité, je souhaite utiliser la diffusion numpy plutôt qu'une boucle Python.

Pour un exemple reproductible minimum, imaginez les trois quantiles que je veux évaluer sont [3.0, 2.0, 1.0], les valeurs moyennes correspondantes sont non tronquées [6.0, 5.0, 4.0], la coupure est inférieure à 1.5, et l'écart-type est non tronquée 3.0.

Les évaluer individuellement fonctionne comme prévu. Si je lance

import numpy as np 
from scipy.stats import truncnorm 
print truncnorm.logpdf(3.0, a=(1.5-6.0)/3.0, b=np.inf, loc=6.0, scale=3.0) 
print truncnorm.logpdf(2.0, a=(1.5-5.0)/3.0, b=np.inf, loc=5.0, scale=3.0) 
print truncnorm.logpdf(1.0, a=(1.5-4.0)/3.0, b=np.inf, loc=4.0, scale=3.0) 

Je reçois

-2.44840736626 
-2.3878150686 
-inf 

(La dernière valeur est -inf parce 1.0 est inférieure à la coupure). L'utilisation de la diffusion numpy pour deux valeurs à la fois fonctionne également comme prévu. Si je lance

print truncnorm.logpdf(
    np.array([3.0, 2.0]), 
    a=(1.5-np.array([6.0, 5.0]))/3.0, 
    b=np.inf, 
    loc=np.array([6.0, 5.0]), 
    scale=3.0 
) 
print truncnorm.logpdf(
    np.array([2.0, 1.0]), 
    a=(1.5-np.array([5.0, 4.0]))/3.0, 
    b=np.inf, 
    loc=np.array([5.0, 4.0]), 
    scale=3.0 
) 

Je reçois

[-2.44840737 -2.38781507] 
[-2.38781507  -inf] 

Cependant, si je tente d'évaluer trois valeurs à la fois en exécutant:

print truncnorm.logpdf(
    np.array([3.0, 2.0, 1.0]), 
    a=(1.5-np.array([6.0, 5.0, 4.0]))/3.0, 
    b=np.inf, 
    loc=np.array([6.0, 5.0, 4.0]), 
    scale=3.0 
) 

Je reçois une erreur:

Traceback (most recent call last): 
    File "truncnorm_error.py", line 25, in <module> 
    scale=3.0 
    File "C:\Python27\lib\site-packages\scipy\stats\_distn_infrastructure.py", line 1701, in logpdf 
    place(output, cond, self._logpdf(*goodargs) - log(scale)) 
    File "C:\Python27\lib\site-packages\scipy\stats\_continuous_distns.py", line 4853, in _logpdf 
    return _norm_logpdf(x) - self._logdelta 
ValueError: operands could not be broadcast together with shapes (2,) (3,) 

Qu'est-ce qui me manque? J'utilise Python 2.7, numpy 1.13, et scipy 0.19.

Answer 1

La raison pour laquelle cela ne fonctionne pas, car logpdf vérifie les quantiles pour s'assurer qu'ils sont plus grands que le seuil. Si vous avez une valeur inférieure à la troncature, apparemment, cela fonctionne pour les tailles 1 et 2, mais pas pour 3. Donc, cela peut être le bug.

Si vous fournissez des valeurs supérieures à la troncature, cela fonctionne correctement. Par exemple, cela fonctionne, où j'ai changé 1.0 à 1.6 pour les quantiles:

print truncnorm.logpdf(
    np.array([3.0, 2.0, 1.6]), 
    a=(1.5-np.array([6.0, 5.0, 4.0]))/3.0, 
    b=np.inf, 
    loc=np.array([6.0, 5.0, 4.0]), 
    scale=3.0) 

Answer 2

Merci, tous. En attendant, je roulais mon:

def left_truncnorm_logpdf(x, untruncated_mean, untruncated_std_dev, left_cutoff): 
    f = np.array(np.subtract(stats.norm.logpdf(x, loc=untruncated_mean, scale=untruncated_std_dev), 
          np.log(1 - stats.norm.cdf(left_cutoff, loc=untruncated_mean, scale=untruncated_std_dev)))) 
    f[x < left_cutoff] = -np.inf 
    return f 

Il est inélégant et je suis sûr qu'il a des problèmes, mais il semble fonctionner pour mes besoins (par exemple, il diffuse correctement des arguments de vecteur pour x et untruncated_mean).