R语言 如何写分母求和?

piwo6bdm  于 2022-12-06  发布在  其他
关注(0)|答案(1)|浏览(151)

我首先对一个随机矩阵及其特征值和特征向量进行采样。

n <- 2000

#Sample GOE random matrix
A <- matrix(rnorm(n*n, mean=0, sd=1), n, n) 
G <- (A + t(A))/sqrt(2*n)
ev <- eigen(G)
l <- ev$values
v <- ev$vectors

我还生成了均匀分布在球面上的随机向量

#size of multivariate distribution
mean <- rep(0, n) 
var <- diag(n)

#make this example reproducible
set.seed(101)

#simulate bivariate normal distribution
initial <- mvrnorm(n=10, mu=mean, Sigma=var)

#normalized the first possible initial value, the initial data uniformly distributed on the sphere
x_0 <- initial[1, ]/norm(initial[1, ], type="2")

所以x_0是一个均匀分布在球面上的随机向量,v的每一列都是矩阵的特征向量。
现在我想定义以下函数

我的代码如下。但我的问题是如何写分母。

h10 <- x_0 %*% v[, n]

sum((x_0 %*% v[, i])^2) #How to fix this one?

h1 <- function(t) {
  abs(h10)*exp(-2*l[n]*t)/(sqrt())
}
k3bvogb1

k3bvogb11#

这里有一个方法,把计算分解成numerdenom,利用R的矢量化能力来计算矩阵乘积、幂和平方根。
这个函数是t上的标量函数,所以我也对"t"上的标量函数进行了矢量化(与上面的无关)。

#n <- 2000
n <- 20

#Sample GOE random matrix
A <- matrix(rnorm(n*n, mean=0, sd=1), n, n) 
G <- (A + t(A))/sqrt(2*n)
ev <- eigen(G)
l <- ev$values
v <- ev$vectors

#size of multivariate distribution
mean <- rep(0, n) 
var <- diag(n)

#make this example reproducible
set.seed(101)

#simulate bivariate normal distribution
initial <- MASS::mvrnorm(n=10, mu=mean, Sigma=var)

#normalized the first possible initial value, the initial data uniformly distributed on the sphere
x_0 <- initial[1, ]/norm(initial[1, ], type="2")

h1_scalar <- function(t, x0, v, lambda) {
  h10 <- x0 %*% v[, n]
  numer <- abs(h10) * exp(-2*lambda[n] * t)
  h <- sum((x0 %*% v)^2)
  denom <- h * exp(-4*lambda * t)
  c(numer)/sqrt(denom)
}
h1 <- Vectorize(h1_scalar, "t")

h1(1, x_0, v, l)
#>              [,1]
#>  [1,] 13.62478196
#>  [2,] 12.53497558
#>  [3,]  7.50725898
#>  [4,]  5.66588613
#>  [5,]  4.22446682
#>  [6,]  2.71595400
#>  [7,]  1.97947406
#>  [8,]  1.13463313
#>  [9,]  0.91133245
#> [10,]  0.75639112
#> [11,]  0.57187518
#> [12,]  0.36569767
#> [13,]  0.29825876
#> [14,]  0.20393953
#> [15,]  0.18553709
#> [16,]  0.12584079
#> [17,]  0.08098017
#> [18,]  0.04032765
#> [19,]  0.03189976
#> [20,]  0.01920207

h1(1:10, x_0, v, l)
#>              [,1]         [,2]         [,3]         [,4]         [,5]
#>  [1,] 13.62478196 9.667431e+03 6.859502e+06 4.867143e+09 3.453470e+12
#>  [2,] 12.53497558 8.182744e+03 5.341637e+06 3.486983e+09 2.276278e+12
#>  [3,]  7.50725898 2.935045e+03 1.147488e+06 4.486229e+08 1.753940e+11
#>  [4,]  5.66588613 1.671813e+03 4.932958e+05 1.455550e+08 4.294841e+10
#>  [5,]  4.22446682 9.293852e+02 2.044653e+05 4.498249e+07 9.896175e+09
#>  [6,]  2.71595400 3.841464e+02 5.433393e+04 7.685029e+06 1.086976e+09
#>  [7,]  1.97947406 2.040570e+02 2.103553e+04 2.168479e+06 2.235409e+08
#>  [8,]  1.13463313 6.704446e+01 3.961597e+03 2.340872e+05 1.383201e+07
#>  [9,]  0.91133245 4.325194e+01 2.052742e+03 9.742338e+04 4.623725e+06
#> [10,]  0.75639112 2.979510e+01 1.173662e+03 4.623188e+04 1.821126e+06
#> [11,]  0.57187518 1.703156e+01 5.072332e+02 1.510640e+04 4.498980e+05
#> [12,]  0.36569767 6.964603e+00 1.326388e+02 2.526066e+03 4.810817e+04
#> [13,]  0.29825876 4.632745e+00 7.195874e+01 1.117709e+03 1.736097e+04
#> [14,]  0.20393953 2.165982e+00 2.300425e+01 2.443214e+02 2.594866e+03
#> [15,]  0.18553709 1.792724e+00 1.732192e+01 1.673705e+02 1.617192e+03
#> [16,]  0.12584079 8.246978e-01 5.404658e+00 3.541943e+01 2.321213e+02
#> [17,]  0.08098017 3.415146e-01 1.440257e+00 6.073940e+00 2.561540e+01
#> [18,]  0.04032765 8.469502e-02 1.778741e-01 3.735663e-01 7.845536e-01
#> [19,]  0.03189976 5.299400e-02 8.803717e-02 1.462532e-01 2.429656e-01
#> [20,]  0.01920207 1.920207e-02 1.920207e-02 1.920207e-02 1.920207e-02
#>               [,6]         [,7]         [,8]         [,9]        [,10]
#>  [1,] 2.450401e+15 1.738676e+18 1.233673e+21 8.753500e+23 6.211024e+26
#>  [2,] 1.485938e+15 9.700100e+17 6.332157e+20 4.133587e+23 2.698377e+26
#>  [3,] 6.857222e+13 2.680906e+16 1.048129e+19 4.097776e+21 1.602070e+24
#>  [4,] 1.267263e+13 3.739268e+15 1.103333e+18 3.255564e+20 9.606077e+22
#>  [5,] 2.177164e+12 4.789775e+14 1.053753e+17 2.318264e+19 5.100194e+21
#>  [6,] 1.537426e+11 2.174546e+13 3.075693e+15 4.350281e+17 6.153068e+19
#>  [7,] 2.304404e+10 2.375530e+12 2.448850e+14 2.524434e+16 2.602350e+18
#>  [8,] 8.173208e+08 4.829476e+10 2.853694e+12 1.686222e+14 9.963737e+15
#>  [9,] 2.194425e+08 1.041477e+10 4.942861e+11 2.345887e+13 1.113361e+15
#> [10,] 7.173619e+07 2.825769e+09 1.113102e+11 4.384635e+12 1.727157e+14
#> [11,] 1.339884e+07 3.990437e+08 1.188430e+10 3.539378e+11 1.054096e+13
#> [12,] 9.162056e+05 1.744886e+07 3.323084e+08 6.328713e+09 1.205284e+11
#> [13,] 2.696616e+05 4.188555e+06 6.505930e+07 1.010542e+09 1.569639e+10
#> [14,] 2.755930e+04 2.926992e+05 3.108672e+06 3.301629e+07 3.506563e+08
#> [15,] 1.562587e+04 1.509826e+05 1.458846e+06 1.409588e+07 1.361993e+08
#> [16,] 1.521207e+03 9.969234e+03 6.533339e+04 4.281625e+05 2.805963e+06
#> [17,] 1.080269e+02 4.555776e+02 1.921290e+03 8.102585e+03 3.417073e+04
#> [18,] 1.647698e+00 3.460449e+00 7.267539e+00 1.526308e+01 3.205510e+01
#> [19,] 4.036306e-01 6.705381e-01 1.113943e+00 1.850556e+00 3.074266e+00
#> [20,] 1.920207e-02 1.920207e-02 1.920207e-02 1.920207e-02 1.920207e-02

创建于2022年11月28日,使用reprex v2.0.2

相关问题