Python中用于大量观察的多元t分布的密度答案

【问题标题】：Density of multivariate t distribution in Python for large number of observationsPython中用于大量观察的多元t分布的密度
【发布时间】：2017-04-22 05:47:00
【问题描述】：

我正在尝试评估 13 维向量的多元 t 分布的密度。使用 R 中 mvtnorm 包中的 dmvt 函数，我得到的结果是

[1] 1.009831e-13

当我尝试自己用 Python 编写函数时（感谢这篇文章中的建议： multivariate student t-distribution with python)，我意识到 gamma 函数的值非常高（鉴于我有 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 等效函数相同的结果吗？

使用 as x = (0,0) 会产生类似的结果（在一定程度上，舍弃四舍五入），但使用 x = (1,1)1 我会得到显着的不同！

【问题讨论】：

尝试在一些不太极端的输入上运行你的函数，例如二维的 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 结果是 0.03062938 和0.0530516476973 在 Python 中

标签： python r python-2.7 statistics

【解决方案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))

【讨论】：