Pythonの多変量分布の観測多数の観測の場合

Question

私は13-dベクトルの多変量t分布の密度を評価しようとしています。 、私は気づい：Rにおけるmvtnormパッケージからdmvt機能を使用して、私が得る結果は、私はPythonで自分で（ multivariate student t-distribution with pythonこの記事での提案のおかげで）関数を記述しようとした

[1] 1.009831e-13 

がありますガンマ関数は非常に高い値をとりました（私はn = 7512の観測値を持っています）、関数が範囲外になりました。

私は対数スケールにそれを変換するためにmath.lgamma（）とnp.linalg.slogdet（）関数を使用して、アルゴリズムを変更しようとしたが、私が得た結果はこれがある

8.97669876e-15 

ました私がPythonで使用した関数は以下の通りです：

def dmvt(x,mu,Sigma,df,d): 
    ''' 
    Multivariate t-student density: 
    output: 
     the density of the given element 
    input: 
     x = parameter (d dimensional numpy array or scalar) 
     mu = mean (d dimensional numpy array or scalar) 
     Sigma = scale matrix (dxd numpy array) 
     df = degrees of freedom 
     d: dimension 
    ''' 
    Num = math.lgamma(1. *(d+df)/2) - math.lgamma(1.*df/2) 
    (sign, logdet) = np.linalg.slogdet(Sigma) 
    Denom =1/2*logdet + d/2*(np.log(pi)+np.log(df)) + 1.*((d+df)/2)*np.log(1 + (1./df)*np.dot(np.dot((x - mu),np.linalg.inv(Sigma)), (x - mu))) 
    d = 1. * (Num - Denom) 
    return np.exp(d) 

この関数がRと同等の結果を出さない理由は何ですか？

x = (0,0)と同様の結果が得られます（一点まで丸めます）が、x = (1,1) 1と大きな違いがあります。

Answer 1

最後に、Rのmvtnormパッケージのコードを翻訳して管理しました。次のスクリプトは数値アンダーフローなしで動作します。

import numpy as np 
import scipy.stats 
import math 
from math import lgamma 
from numpy import matrix 
from numpy import linalg 
from numpy.linalg import slogdet 
import scipy.special 
from scipy.special import gammaln 

mu = np.array([3,3]) 
x = np.array([1, 1]) 
Sigma = np.array([[1, 0], [0, 1]]) 
p=2 
df=1 

def dmvt(x, mu, Sigma, df, log): 
    ''' 
    Multivariate t-student density. Returns the density 
    of the function at points specified by x. 

    input: 
     x = parameter (n x d numpy array) 
     mu = mean (d dimensional numpy array) 
     Sigma = scale matrix (d x d numpy array) 
     df = degrees of freedom 
     log = log scale or not 

    ''' 
    p = Sigma.shape[0] # Dimensionality 
    dec = np.linalg.cholesky(Sigma) 
    R_x_m = np.linalg.solve(dec,np.matrix.transpose(x)-mu) 
    rss = np.power(R_x_m,2).sum(axis=0) 
    logretval = lgamma(1.0*(p + df)/2) - (lgamma(1.0*df/2) + np.sum(np.log(dec.diagonal())) \ 
     + p/2 * np.log(math.pi * df)) - 0.5 * (df + p) * math.log1p((rss/df)) 
    if log == False:  
     return(np.exp(logretval)) 
    else: 
     return(logretval) 


print(dmvt(x,mu,Sigma,df,True)) 
print(dmvt(x,mu,Sigma,df,False))