acos(1) возвращает NaN для некоторых значений, а не для других

Question

У меня есть список значений широты и долготы, и я пытаюсь найти расстояние между ними. Используя стандартный метод большого круга, мне нужно найти:

acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1))

И умножьте это на радиус Земли, в единицах, которые я использую. Это действительно до тех пор, пока значения, которые мы принимаем, находятся в диапазоне [-1,1]. Если они даже немного выходят за пределы этого диапазона, он вернется NaNдаже если разница связана с округлением.

У меня проблема в том, что иногда, когда два значения lat/long совпадают, это дает мне NaN ошибка. Не всегда, даже для одной и той же пары чисел, но всегда одни и те же в списке. Например, я остановил человека на дороге в пустыне:

Time  |lat     |long
1:00PM|35.08646|-117.5023
1:01PM|35.08646|-117.5023
1:02PM|35.08646|-117.5023
1:03PM|35.08646|-117.5023
1:04PM|35.08646|-117.5023

Когда я вычисляю расстояние между последовательными точками, третье значение, например, всегда будет NaNхотя другие нет. Это кажется странной ошибкой с округлением R.

Answer 1

Не могу точно сказать, не видя ваших данных (попробуйте dput), но это, скорее всего, следствие FAQ 7.31.

(x1 <- 1)
## [1] 1
(x2 <- 1+1e-16)
## [1] 1
(x3 <- 1+1e-8)
## [1] 1
acos(x1)
## [1] 0
acos(x2)
## [1] 0
acos(x3)
## [1] NaN

То есть, даже если ваши значения настолько похожи, что их печатные представления одинаковы, они все равно могут отличаться: некоторые будут в пределах .Machine$double.eps а другие не будут...

Один из способов убедиться, что входные значения ограничены [-1,1], это использовать pmax а также pmin: acos(pmin(pmax(x,-1.0),1.0))

Answer 2

Простой обходной путь - использовать pmin(), например так:

acos(pmin(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1),1))

Это теперь гарантирует, что потеря точности приводит к значению не выше, чем ровно 1.

Это не объясняет, что происходит, однако.

(Правка: Мэтью Лундберг указал, что мне нужно использовать pmin для работы с векторизованными входами. Это решает проблему с получением его работы, но я все еще не уверен, почему он округляется неправильно.)

Answer 3

Я только что столкнулся с этим. Это вызвано тем, что ввод больше 1. Из-за вычислительной ошибки мой внутренний продукт между единичными нормами становится немного больше 1 (например, 1 + 0,00001). А также acos()может иметь дело только с [-1,1]. Таким образом, мы можем ограничить верхнюю границу ровно 1, чтобы решить проблему.

Для numpy: np.clip(your_input, -1, 1)

Для Питорча: torch.clamp(your_input, -1, 1)