R:矩阵乘法错误(不一致的参数)

Question

我正在使用 R 编程语言。

我有以下数据:

1) 均值向量(4 行，1 列)

4 个变量(x1、x2、x3、x4)

      5.0060022 
     3.4280049 
    1.4620007 
   0.2459998

2)协方差矩阵(4行4列)

4 个变量(对角线元素是 x1、x2、x3、x4，成对元素是例如第二个元素:(x1,x2)、第三个元素 (x1,x3)、第四个元素 (x1, x4) 等)

   0.15065114  0.13080115   0.02084463  0.01309107
   0.13080115  0.17604529   0.01603245  0.01221458
  0.02084463  0.01603245   0.02808260  0.00601568
  0.01309107  0.01221458   0.00601568  0.01042365

问题:我想获取上述数据并创建一个函数(有 4 个输入:x1、x2、x3、x4 和一个数字作为输出)，格式如下:

这是我到目前为止尝试过的:

my_function <- function(x_one, x_two, x_three, x_four)

{


sigma1.pre <- c(0.15065114 , 0.13080115 ,  0.02084463 , 0.01309107 , 0.13080115 , 0.17604529 ,  0.01603245 , 0.01221458 , 0.02084463 , 0.01603245  , 0.02808260 , 0.00601568 , 0.01309107 , 0.01221458 ,  0.00601568 , 0.01042365)
sigma1 <- matrix(sigma1.pre, nrow=4, ncol= 4, byrow = TRUE)
sigma1_inv <- ginv(sigma1)
det_sigma1_inv <- det(sigma1_inv)
denom = sqrt( (2*pi)^4 * det_sigma1_inv) 

x_one = x1 - 5
x_two = x2 - 3.42
x_three = x3 - 1.462
x_four = x4 - 0.245

x_t = c(x_one, x_two, x_three, x_four)
x_t_one <- matrix(x_t, nrow=4, ncol= 1, byrow = TRUE)
x_t_one_t = -0.5 * t(x_t_one)
x_t_two =  matrix(x_t, nrow=1, ncol= 4, byrow = TRUE)

num = exp(x_t_two  %*%  sigma1_inv  %*%  x_t_one_t)

answer = num/denom

return(answer)
}

问题:当我尝试运行这个函数时:

my_function(1,2,3,4)

我收到以下错误:

Error in x_t_two %*% sigma1_inv %*% x_t_one_t : non-conformable arguments

我认为错误是由于矩阵乘法而发生的

num = exp(x_t_two  %*%  sigma1_inv  %*%  x_t_one_t)

我试图改变矩阵乘法的顺序:

num = exp( x_t_one_t   %*%  sigma1_inv  %*% x_t_two )

但是错误依旧存在。

有人可以告诉我如何解决这个问题吗？

谢谢!

引用资料:

• https://stat.ethz.ch/R-manual/R-devel/library/base/html/matmult.html
• https://en.wikipedia.org/wiki/Multivariate_normal_distribution
• https://wikimedia.org/api/rest_v1/media/math/render/svg/c66e6f6abd66698181e114a4b00da97446efd3c4

Answer 1

正如我上面提到的，dmvnorm 函数返回您显示的函数的值。

dmvnorm(c(5,3,1,0),m,v)
[1] 0.01074766

这是我的手动版本，

func <- function(vec, m, v){
  if (length(vec) != length(m)) {
    stop("dimension error")
  } # and several more
  a <- t(vec - m) %*% solve(v) %*% (vec - m)
  k <- length(vec)
  return(exp(-a/2)/sqrt((2*pi)^k * det(v)))
}

func(c(5,3,1,0) , m, v)
           [,1]
[1,] 0.01074766

在您的函数中，您的函数不起作用的主要原因是行 num = exp(x_t_two %*% sigma1_inv %*% x_t_one_t)，x_t_one_t 是错误的。当您将其设置为 nrow = 4, ncol = 1 时，它已经是 4*1，您不需要转置它。我对您的功能添加了更多评论。

 my_function <- function(x_one, x_two, x_three, x_four)
  
{
  
  
  sigma1.pre <- c(0.15065114 , 0.13080115 ,  0.02084463 , 0.01309107 , 0.13080115 , 0.17604529 ,  0.01603245 , 0.01221458 , 0.02084463 , 0.01603245  , 0.02808260 , 0.00601568 , 0.01309107 , 0.01221458 ,  0.00601568 , 0.01042365)
  sigma1 <- matrix(sigma1.pre, nrow=4, ncol= 4, byrow = TRUE)
  # You can also use solve instead of ginv, solve is in base R
  sigma1_inv <- ginv(sigma1)
  det_sigma1_inv <- det(sigma1_inv)
  # In here, not det_sigma1_inv, just use det(sigma1) will work.
  denom = sqrt( (2*pi)^4 * det(sigma1)) 
  
  #in below part, I recommend another way.
  #m <- c( 5.0060022, 3.4280049, 1.4620007, 0.2459998)
  #x_t = c(x_one, x_two, x_three, x_four)
  #There was no input x1, x2, x3, x4
  x_one = x_one - 5.0060022
  x_two = x_two - 3.4280049
  x_three = x_three - 1.4620007
  x_four = x_four - 0.2459998
  
  
  # Vectors and matrices are handle as vector and matrices. You do not need to 
  #change vectors to matrices.
  #x_t_t = x_t - m
  x_t = c(x_one, x_two, x_three, x_four)
  x_t_one <- matrix(x_t, nrow=4, ncol= 1, byrow = TRUE)
  x_t_two =  matrix(x_t, nrow=1, ncol= 4, byrow = TRUE)
  
  
  # In this part, as it's (x-mu)^T * SIGMA * (x-mu), dimension of x_t_one_t was wrong
  # You may try another way.
  #num = exp(-0.5 * t(x_t_t) %*% sigma1_inv %*% x_t_t)
  num = exp(-0.5 * x_t_two  %*%  sigma1_inv  %*%  x_t)
  
  
  answer = num/denom
  
  return(answer)
}
my_function(5,3,1,0)
           [,1]
[1,] 0.01074766