Java实现多元t分布函数(Multivariate t distributions)

短命女 2022-01-22 08:27 343阅读 0赞

本文作者：合肥工业大学 管理学院 钱洋 email：1563178220@qq.com 内容可能有不到之处，欢迎交流。  
**未经本人允许禁止转载**

### 文章目录 ###

*  问题背景
 *  编程实现
 *  使用方法

# 问题背景 #

最近，在独立实现Gaussian LDA算法时，遇到了Multivariate t distributions。Gaussian LDA对应的论文是：  
Das R, Zaheer M, Dyer C. Gaussian lda for topic models with word embeddings\[C\]//Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers). 2015, 1: 795-804.

我的另外一篇博客，已对参数学习公式的推理进行了详细说明，即：

https://qianyang-hfut.blog.csdn.net/article/details/79528781

其中，在实现该算法时，需要对每个单词所属的主题的概率进行计算，即通过如下公式：  
![在这里插入图片描述][2019060517174274.png_pic_center]

从上面的公式可以看到，最为主要的是计算t分布函数，其输入是word2Vec产生的词向量，以及均值向量和协方差矩阵。

**多元t分布函数的公式如下**：  
![在这里插入图片描述][20190605171947558.png_pic_center]

# 编程实现 #

在math3中，提供了很多函数，如Gamma函数，高斯函数等；也提供了很多分布，如高斯分布，Beta分布，二项分布等的采样。下面使用math3实现一元t分布的采样：

import org.apache.commons.math3.distribution.TDistribution;
    
    public class Test3 { 
    
    	public static void main(String[] args) { 
    		//定义自由度
    		TDistribution t = new TDistribution(5);
    		for (int i = 0; i < 100; i++) { 
    			System.out.println(t.sample());
    		}
    	}
    }

程序的输入结果为：

![在这里插入图片描述][2019060517284718.png_pic_center]

下图为一元t分布对应的分布图。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9xaWFueWFuZy1oZnV0LmJsb2cuY3Nkbi5uZXQ_size_16_color_FFFFFF_t_70_pic_center]  
**在math3只提供了一元t分布和其概率密度函数，但并多元t分布的函数。**

在实现Gaussian LDA的过程中，我自己又编写了多元t分布的函数，即下面公式对应的代码：  
![在这里插入图片描述][20190605173230162.png_pic_center]

如下为Java代码：

/**The probability density function of the Multivariate t distributions * Reference: Conjugate Bayesian analysis of the Gaussian distribution * * * @param ArrayRealVector dataPoint * @param ArrayRealVector meansVector * @param RealMatrix covarianceMatrix * @param double degreesOfFreedom * @return The log probability value * * * @author Qianyang * ****/
    	public static double logMultivariateTDensity(ArrayRealVector dataPoint, ArrayRealVector meansVector, RealMatrix covarianceMatrix, double degreesOfFreedom){ 
    		LUDecomposition covariance = new LUDecomposition(covarianceMatrix);
    		double logprob_left = Gamma.logGamma((degreesOfFreedom + dataPoint.getDimension())/2.0) - 
    				(Gamma.logGamma(degreesOfFreedom / 2.0) + 0.5 * Math.log(covariance.getDeterminant()) + 
    						dataPoint.getDimension()/2.0 * (Math.log(degreesOfFreedom) + Math.log(Math.PI)));		
    		// compute x-u
    		ArrayRealVector var = dataPoint.add(meansVector.mapMultiplyToSelf(-1.0));
    		// (x-u) to matrix
    		RealMatrix realMatrix = new Array2DRowRealMatrix(var.getDataRef());
    		//compute left
    		double logprob_right = Math.log(1 + realMatrix.transpose().multiply(new LUDecomposition(covarianceMatrix).getSolver().getInverse())
    				.multiply(realMatrix).getData()[0][0]/degreesOfFreedom);
    		return Math.exp(logprob_left -(degreesOfFreedom + dataPoint.getDimension())/2.0 * logprob_right);
    	}

# 使用方法 #

如下为使用方法：

ArrayRealVector dataPoint = new ArrayRealVector(new double[]{ 2,5});
    		ArrayRealVector meansVector = new ArrayRealVector(new double[]{ 2,1});
    		RealMatrix covarianceMatrix =  MatrixUtils.createRealMatrix(new double[][]{ { 1, 0}, { 0, 1}});
    		double logt = logMultivariateTDensity(dataPoint,meansVector,covarianceMatrix,10);
    		System.out.println(logt);

主方法执行这些代码，会输出如下结果。  
![在这里插入图片描述][20190606084546159.png_pic_center]

改写协方差和自由度：

ArrayRealVector dataPoint = new ArrayRealVector(new double[]{ 2,5});
    		ArrayRealVector meansVector = new ArrayRealVector(new double[]{ 2,1});
    		RealMatrix covarianceMatrix =  MatrixUtils.createRealMatrix(new double[][]{ { 3, 0.8}, { 0.8, 5}});
    		double logt = logMultivariateTDensity(dataPoint,meansVector,covarianceMatrix,17);
    		System.out.println(logt);

程序输出结果为：  
![在这里插入图片描述][20190606085018404.png_pic_center]

[2019060517174274.png_pic_center]: /images/20220122/64a6c3a3b10c4892af1bc7dffaa1eac1.png
[20190605171947558.png_pic_center]: /images/20220122/a3d66671689d48819af6f113b9112383.png
[2019060517284718.png_pic_center]: /images/20220122/59c516fe3d2047f7a831e96c5c37b12f.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9xaWFueWFuZy1oZnV0LmJsb2cuY3Nkbi5uZXQ_size_16_color_FFFFFF_t_70_pic_center]: /images/20220122/d3abebb22ebe419ba7ccee03cbcead41.png
[20190605173230162.png_pic_center]: /images/20220122/be2738fbe67a4a52a8c33da8c8538f37.png
[20190606084546159.png_pic_center]: /images/20220122/0cbffb1676944cdc8c2e8dfa18709914.png
[20190606085018404.png_pic_center]: /images/20220122/e3806c58dc8f49feadfbfeadd6d25c89.png