## Background
We're given a coin that shows heads with an unknown probability, $\lambda$. The goal is to use that coin (and possibly also a fair coin) to build a "new" coin that shows heads with a probability that depends on $\lambda$, call it $f(\lambda)$. This is the _Bernoulli factory problem_, and it can be solved only for certain functions $f$. (For example, flipping the coin twice and taking heads only if exactly one coin shows heads, we can simulate the probability $2\lambda(1-\lambda)$.)
This question is a continuation of [another question of mine](https://mathoverflow.net/questions/379858/what-are-ways-to-compute-polynomials-that-converge-from-above-and-below-to-a-con), which seeks ways to compute polynomials that converge from above and below to $f$ in a manner that solves the Bernoulli factory problem for $f$. These polynomials form an _approximation scheme_ for $f$. See that question for more details, but here is a recap:
**Formal statement:** There exist two sequences of polynomials, namely—
- $g_{n}(\lambda): =\sum_{k=0}^{n}a(n, k){n \choose k}\lambda^{k}(1-\lambda)^{n-k}$, and
- $h_{n}(\lambda): =\sum_{k=0}^{n}b(n, k){n \choose k}\lambda^{k}(1-\lambda)^{n-k}$,
for every integer $n\ge1$, such that—
1. $0\le a(n, k)\le b(n, k)\le1$,
2. $\lim_{n}g_{n}(\lambda)=\lim_{n}h_{n}(\lambda)=f(\lambda)$ for every $\lambda\in[0,1]$, and
3. for every $m<n$, both $(g_{n} - g_{m})$ and $(h_{m} - h_{n})$ have non-negative coefficients once $g_{n}$, $g_{m}$, $h_{n}$, and $h_{m}$ are rewritten as degree-$n$ polynomials in Bernstein form,
where $f(\lambda)$ is continuous on $[0, 1]$ (Nacu and Peres 2005; Holtz et al. 2011), and the goal is to find the appropriate values for $a(n, k)$ and $b(n, k)$.
It is allowed for $a(n, k)\lt0$ for a given $n$ and some $k$, in which case all $a(n, k)$ for that $n$ are taken to be 0 instead. It is allowed for $b(n, k)\gt1$ for a given $n$ and some $k$, in which case all $b(n, k)$ for that $n$ are taken to be 1 instead.
## New coins from old, smoothly
This time, though, we focus on a **specific approximation scheme**, the one **presented by [Holtz et al. 2011](https://link.springer.com/content/pdf/10.1007/s00365-010-9108-5.pdf), in the paper "New coins from old, smoothly"**.
The scheme involves building polynomials that are shifted upward and downward to approximate $f$ from above and below.
The scheme achieves a convergence rate that generally depends on the smoothness of $f$; in fact, it can achieve the highest convergence rate possible for functions with that smoothness. Roughly speaking, the rate is $O((1/n)^{2(r+\alpha)})$ when $\lambda$ is close to 0 or 1, and $O((1/n)^{r+\alpha})$ elsewhere.
The scheme is as follows:
Let $f$ be a continuous and $r$-times differentiable function—
- that maps [0, 1] to the open interval (0, 1), and
- whose $r$th derivative is $\alpha$-Hölder continuous.
Let $b = 2^s$ and let $s\ge0$ be an integer. Let $Q_{n, r}f$ be a degree $n+r$ approximating polynomial called a _Lorentz operator_ (see the paper for details on the Lorentz operator). Let $n_0$ be the smallest $n$ such that $Q_{n_0, r}f$ has coefficients within [0, 1]. Define the following for every integer $n \ge n_0$ divisible by $n_{0}b$:
- $f_{n_0} = Q_{n_0, r}f$.
- $f_{n} = f_{n/b} + Q_{n, r}(f-f_{n/b})$ for each integer $n > n_0$.
- $\phi(n, \alpha, \lambda) = \frac{\theta_{\alpha}}{n^{\alpha}}+(\frac{\lambda(1-\lambda)}{n})^{\alpha/2}$.
Then the coefficients for the degree $n+r$ polynomial that approximates $f$ are—
- $g(n, r, k) = f_{n}(k/(n+r)) - D * B_{n, \phi(n, r+\alpha, \lambda)}(k/(n+r))$, and
- $h(n, r, k) = f_{n}(k/(n+r)) + D * B_{n, \phi(n, r+\alpha, \lambda)}(k/(n+r))$,
where $B_{n, F}$ is the degree-$n$ Bernstein polynomial of $F$, which I know how to compute.
The problem here is:
- This method uses three constants, namely $s$, $\theta_{\alpha}$, and $D$, that are vaguely defined in the paper. For example, the paper says only that $D$ should be chosen "large enough".
- The paper has no examples of how the scheme works for a selection of functions $f$, making it hard to understand and apply the Holtz method (especially in an automated way).
## Questions
1. What are practical lower bounds for $s$, $\theta_{\alpha}$, and $D$ for the method above, given a function $f$, with or without additional assumptions on $f$ (such as being $C^2$ continuous, concave, convex, piecewise polynomial, real analytic, or any combination of these)?
2. What are other practical ways to apply the Holtz method to certain functions to find polynomials that converge to those functions?
3. To help me and others (especially programmers) understand and apply the Holtz method, can you provide examples of the method (in terms of Bernstein polynomial approximation schemes) for various functions or classes of functions (such as $C^2$ continuous, concave, convex, piecewise polynomial, real analytic, or any combination of these)?
## My Attempt
The following Python code shows my attempt to implement the Holtz approximation scheme using the SymPy computer algebra library. It includes an implementation of the Lorentz operator as well as a method to find Bernstein polynomials for a function.
```
TAURJ_CACHE = {}
TNJ_CACHE = {}
TAURJ_SYMBOL = symbols("T")
def taurj_x(r, j, n):
if j == 0:
return S(1)
if j == 1:
return S(0)
if (r, j, n) in TAURJ_CACHE:
return TAURJ_CACHE[r, j, n]
ret = 0
for l in range(2, j + 1):
tnj = None
if (n, l) in TNJ_CACHE:
tnj = TNJ_CACHE[(n, l)]
else:
tnj = sum(
(k - n * TAURJ_SYMBOL) ** l
* binomial(n, k)
* TAURJ_SYMBOL ** k
* (1 - TAURJ_SYMBOL) ** (n - k)
for k in range(0, n + 1)
)
ret -= tnj * taurj_x(r - l, j - l, n) / factorial(l)
TAURJ_CACHE[(r, j, n)] = ret.simplify()
return ret
def taurj(r, j, x, n):
if j == 0:
return 1
if j == 1:
return 0
return taurj_x(r, j, n).subs(TAURJ_SYMBOL, x)
def lorentz_n_k(func, x, n, r, k, pt=None):
# kth Bernstein coefficient of the Lorentz operator
# of degree n and smoothness r at the point
# pt (or at x if not given).
if pt == None:
pt = x
ret = 0
for j in range(0, r + 1):
ret += diff(func, (x, j)).subs(x, S(k) / n) * taurj(r, j, pt, n) / n ** j
return ret
def lorentz_n(func, x, n, r, pt=None):
# Lorentz operator.
# Create a polynomial that approximates func, which in turn uses
# the symbol x. The polynomial's degree is n, is assumed to
# be at least r times differentiable, and is evaluated
# at the point pt (or at x if not given).
# NOTE: If r is 0, this is simply the Bernstein operator.
if pt == None:
pt = x
ret = 0
diffs = [diff(func, (x, j)) for j in range(0, r + 1)]
for k in range(0, n + 1):
rm = 0
for j in range(0, r + 1):
rm += diffs[j].subs(x, S(k) / n) * taurj(r, j, pt, n) / n ** j
ret += rm * bernstein_n_k(n, k, pt)
return ret
def bernstein_n_k(n, k, x):
return binomial(n, k) * x ** k * (1 - x) ** (n - k)
```
REFERENCES:
- Keane, M. S., and O'Brien, G. L., "A Bernoulli factory", ACM Transactions on Modeling and Computer Simulation 4(2), 1994.
- Holtz, O., Nazarov, F., Peres, Y., "New Coins from Old, Smoothly", Constructive Approximation 33 (2011).