如何求出时间t>0的值,使得R中的h(t)=0.01?

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

问题中的改进答案:How to get the value of t so that my function h(t)=epsilon for a fixed epsilon?
我的问题是:
考虑一个随机矩阵,对它的特征向量v_i和特征值lambda_i进行采样。给定初始数据x_0,我想得到固定epsilon=0.01的命中时间t_n:=\inf\{t>0: h_1(t)\ge \epsilon\}。这里函数h_1(t)

给出。
我已经为这些设置和函数h_1(t)编写了代码:

#make this example reproducible
set.seed(100001)
n <- 300
#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)

#simulate bivariate normal distribution
initial <- MASS::mvrnorm(n=1000, mu=mean, Sigma=var) #ten random vectors
#normalized the first possible initial value, the initial data uniformly distributed on the sphere
xmats <- lapply(1:1000, function(i) initial[i, ]/norm(initial[i, ], type="2"))

h1t <- function(t,x_0) {
  h10 <- c(x_0 %*% v[, n])
  denom <- vapply(t, function(.t) {
    sum((x_0 %*% v)^2 * exp(-4*(l - l[n]) * .t))
  }, numeric(1L))
  abs(h10) / sqrt(denom)
}

我用这个问题的答案来计算击球时间:

find_t <- function(x, epsilon = 0.01, range = c(-50, 50)) {
  uniroot(function(t) h1t(t, x) - epsilon, range,
          tol = .Machine$double.eps)$root
}

res <- lapply(xmats, find_t)

输出res

[[995]]
[1] -0.2698699

[[996]]
[1] -0.3138642

[[997]]
[1] -0.4417028

[[998]]
[1] -0.04204563

[[999]]
[1] -0.4150783

[[1000]]
[1] -0.3695955

**问题:**但此输出res将包含负值。如何修复此问题?

如果我绘制函数h_1(t)的图形:我们可以看到,对于epsilon=0.01,时间t的值应该是正的......所以,看起来函数find_t有问题。

h1t <- function(t,x_0=unlist(xmats[1000])) {
  h10 <- c(x_0 %*% v[, n])
  denom <- vapply(t, function(.t) {
    sum((x_0 %*% v)^2 * exp(-4*(l - l[n]) * .t))
  }, numeric(1L))
  abs(h10) / sqrt(denom)
}

plot(h1t,0,200)

更新:
我发现,如果我选择n=1000为矩阵的大小,就会出现错误:

Error in uniroot(function(t) h1t(t, x) - epsilon, range, tol = .Machine$double.eps) : 
f() values at end points not of opposite sign
xxls0lw8

xxls0lw81#

你的res没有问题,如下图所示,水平线画在y == epsilon == 0.01上。
你把横坐标当成纵坐标了,就这样。

res <- lapply(xmats, find_t)

curve(h1t, -1, 1, ylim = c(0, 1))
abline(h = 0.01, v = res[[1000]], col = "red", lty = "dashed")

创建于2022年11月29日,使用reprex v2.0.2
严格递增函数定义为t > 0,但

h1t(0)
#> [1] 0.07184164

在它的定义域中,不存在t对其中h1t(t) == 0.01

相关问题