用于大量观察的Python中多变量t分布的密度

我想要评估13维矢量的多变量t分布密度。从mvtnorm包中的R使用dmvt功能，结果我得到的是用于大量观察的Python中多变量t分布的密度

[1] 1.009831e-13

当我试图通过自己在Python写的函数（感谢在这个岗位的建议： multivariate student t-distribution with python），我意识到，伽马函数的值非常高（因为我有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我得到了显着的差异！

来源

2016-12-07 Leonidas Souliotis

尝试在一些不那么极端的输入上运行你的函数，例如什么是2维的dvmt（0）。这会告诉你，你是否正在陷入舍入误差（由于双精度FP算法的限制，8e-15为零），或者如果你的代码有错误。 –

@HongOoi非常感谢您的评论！使用以下输入： 'x =（0,0） mu =（0,0） sigma = diag（2）' 我在R中获得了'0.1591549'，在Python中获得了'0.159154943092'，这似乎是直到一个点但是当我使用x =（1,1）时，R的结果是Python中的“0.03062938”和“0.0530516476973” –

我终于成功地从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))

来源

2017-02-20 09:15:04

用于大量观察的Python中多变量t分布的密度

回答

相关问题