在R中创建大比值 [英] Creating Mills Ratio in R for large values

问题描述

我正在使用R创建一个函数，该函数除其他外还使用Mills Ratio（请参见此处）。这不是一个复杂的公式，起初我只是这样编程：

I'm using R to create a function, that amongst others uses Mills Ratio (See here). This is not a complicated formula, and at first I just programmed it like this:

mill <- function(x) {
  return((1 - pnorm(x)) / dnorm(x))
}

但是我很快发现，对于 x 的非常大的值（x> = 9），该函数返回零。更引人注目的是，它在x> = 37左右开始返回 NaN ，这确实弄乱了我的东西。

I soon found out however, that for very large values (x >= 9) of x , this function returns zero. Even more dramatic, at around x >= 37, it starts returning NaN , which really messes up my stuff.

现在我将函数更改为：

mill <- function(x) {
  if (x >= 9) {
    return(1 / x)
  } else {
    return((1 - pnorm(x)) / dnorm(x))
  }
}

这似乎有效。我使用此函数来计算向量，但是当我使用仿真找到相同的向量时，或多或少都得到相同的答案，只是总是有点偏离。

This seems to work. I use this function to calculate a vector however, and when I use simulation to find the same vector, I get more or less the same answer, only it's always a bit off..

我认为这与我的Mills Ratio的实现有关，因为该函数的其余部分只是指数，R应该没有问题。

I think this has to do with my implementation of Mills Ratio, since the rest of the function is just exponentials, which R should have no trouble with.

我想问你们是否有什么方法可以解决这个问题：要么更好地实现此功能，要么给我另一种找到米尔斯比率的方法（也许通过某种形式的集成，但是我不会遇到同样的问题吗？）。

I want to ask you guys if there is any way to solve this problem: to either implement this function better, or give me another way to find the Mills Ratio (perhaps through integration of some sorts, but wouldn't I run into the same issues there?). Thank you kindly for any help you can provide!

推荐答案

我将对原始的磨机进行两项更改函数。

1. 更改 1-pnorm（x）到 norm（lower.tail = FALSE）


2. 使用 log '

1. Change 1-pnorm(x) to pnorm(lower.tail=FALSE)
2. Use log's and take exponentials if needed.

所以这给出了

new_mill = function(x) 
    pnorm(x, lower.tail=FALSE, log.p=TRUE) - dnorm(x, log=TRUE)

所以

R> exp(new_mill(10))
[1] 0.09903
R> exp(new_mill(40))
[1] 0.02498

将情节视为理智检查

x = seq(0, 10, 0.001)
plot(x, exp(new_mill(x)), type="l")
lines(x, mill(x), col=2)

给予

这篇关于在R中创建大比值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！