Log utility function and the St. Petersburg paradox

Question

In the log utility model the formula solving the St. Petersburg paradox

$$\Delta E(U)=\sum_{k=1}^\infty \frac{1}{2^k}\left [\ln(w + 2^k - c) - \ln(w) \right ]$$

relates the wealth, $w,$ of the individual to how much they should be able to pay to enter the game, $c.$ The $\Delta E(U)$ should be positive.

The question is whether this $c$ value can be calculated in close form (using geometric series), or whether there is a need for numeric optimization.

For instance, the Wikipedia article states

For example, with natural log utility, a millionaire ($\$1,000,000$) should be willing to pay up to $\$20.88,$ a person with $\$1,000$ should pay up to $\$10.95,$ a person with $\$2$ should borrow $\$1.35$ and pay up to $\$3.35.$

How are these values calculated?

---

CODE OF THE FUNCTIONS IN THE ACCEPTED ANSWER IN R STATS:

E = function(w, c, epsilon){
    ans = 0
    k = 1

    while(abs(val <- (log(max(epsilon, w + 2^k - c)) - log(w)) / 2^k) > epsilon){
        k <- k + 1; ans <- ans + val
        if(abs(val) < epsilon) break
    }
    ans
}

find_c = function(w, epsilon=10^(-10)){
    lo = 0
    c = 0
    hi = 10^10

    while(abs(lo - hi) > epsilon){
        c = (hi + lo) / 2
        exp_value = E(w, c, epsilon)

        ifelse(exp_value > 0, lo <- c, hi <- c)
    }
    c
}

find_c(10^6); find_c(10^3)
20.8742701188956
10.9538128679429

Answer 1

My guess would be numerical methods. For example, in Python I get the following results, using a sort of binary search to hone in on the value of $c$ that makes the expected value near-$0$:

from math import log


def E(w, c, epsilon):
    ans = 0
    k = 1

    while True:
        val = (log(max(epsilon, w + 2**k - c)) - log(w)) / 2**k

        if abs(val) < epsilon:
            break

        k += 1
        ans += val

    return ans


def find_c(w, epsilon=10**(-10)):
    lo = 0
    c = 0
    hi = 10**10

    while abs(lo - hi) > epsilon:
        c = (hi + lo) / 2
        exp_value = E(w, c, epsilon)

        if exp_value > 0:
            lo = c
        else:
            hi = c

    return c


for w in [10**6, 10**3, 2]:
    c = find_c(w)

    if w >= c:
        print("For w = {}, c = {:0.2f}".format(w, c))
    else:
        print("For w = {}, borrow {:0.2f} and pay c = {:0.2f}".format(w, c - w, c))

With results:

For w = 1000000, c = 20.87
For w = 1000, c = 10.95
For w = 2, borrow 1.35 and pay c = 3.35