2

In R, I need to calculate sum((v/u) * (d^a))^(1/a) (see below) when d can be from 0 up to 10000 and a may be from 0 to 100. However, this overflows to Inf. v/u is dimensionless, as is a. I therefore assume that the result is on the same scale as d.

Is there a way to avoid the overflow in this calculation, either by some clever math, built-in functions, or by auto-scaling d, and then re-scaling the result without losing precision?

> N <- 200
> a <- 100
> v <- runif(N, min=1, max=20)
> u <- sum(v)
> d <- runif(N, min=1, max=8000)
> sum((v/u) * (d^a))^(1/a)
[1] Inf

enter image description here

M--
  • 25,431
  • 8
  • 61
  • 93
caracal
  • 2,690
  • 19
  • 22
  • 1
    "`a` may be from 0 to 100." and code has `1/a`? To avoid that overflow, treat `a==0` as a special case. Not certain what the return value should be other than some infinity. – chux - Reinstate Monica Nov 14 '17 at 21:36
  • 1
    @chux but here OP defined `a <- 100`. So, that's not the problem. But `d` can be between 1 and 8000 which to the power of 100 can be a very large number. That's the issue I think. Any `d >= 1210` to the power 100 will be `Inf`. – M-- Nov 14 '17 at 22:06
  • https://stackoverflow.com/questions/22466328/how-to-work-with-large-numbers-in-r – M-- Nov 14 '17 at 22:12
  • @chux Right, `a == 0` needs to caught earlier on. The overflow to `Inf` happens for large `a`. – caracal Nov 15 '17 at 09:32

1 Answers1

4

I'm not sure if this is what you look for, but you can rescale your d vector

N <- 200
a <- 100
v <- runif(N, min=1, max=20)
u <- sum(v)
d <- runif(N, min=1, max=8000)
D <- max(d)
D*sum((v/u) * ((d/D)^a))^(1/a) 
[1] 7721.741

Here I used the maximum value, but you can get better results as long as D is large enough to avoid (d/D)^a reaching computer infinity.

The point is, if you choose D too small your sum is ruined because you get Inf in a term, but if you D is chosen excessively large some terms of your sum will be rounded to zero, which is less precise but definitely better than infinity.

VFreguglia
  • 2,129
  • 4
  • 14
  • 35
  • 2
    +1. Since "d" comes from a uniform distribution, using `D=max(d)` is almost perfect. To be completely foolproof, you'd want to catch the bad cases where the v/u term happens to be very small in conjonction with "d" reaching its max.. To prevent that, use this: `D <- max((v/u)^(1/a)*d)` – Fabien Viger Nov 14 '17 at 23:17
  • Thank you Victor and @FabienViger! – caracal Nov 15 '17 at 09:34