Aunque esta es mi pregunta, me voy a dar una respuesta que muestra un posible algoritmo, para conseguir el balanceo de la bola. Una forma de calcular la curva Bézier es De Casteljau del algoritmo, que es numéricamente estable método de cálculo que utiliza un método recursivo relacionados con la propiedad recursiva de la distribución binomial. Este algoritmo puede ser implementado tanto en el estándar de la probabilidad de escala, o en el registro de probabilidad de escala.
Recursiva caracterización de la curva Bézier: Para obtener una relación de recurrencia para la caracterización de esta función, vamos a tomar ventaja de la bien conocida ecuación recursiva:
$$\text{Bin}(x|n, \theta) = (1-\theta) \cdot \text{Bin}(x|n-1, \theta) + \theta \cdot \text{Bin}(x-1|n-1, \theta).$$
El uso de esta ecuación recursiva, para cualquier vector de argumento $\mathbf{f} = (f_0,...,f_n)$, tenemos:
$$\begin{equation} \begin{aligned}
B(\mathbf{f}, \theta)
&= \sum_{x=0}^n f_x \cdot \text{Bin}(x|n,\theta) \\[6pt]
&= \sum_{x=0}^n f_x \Big[ (1-\theta) \cdot \text{Bin}(x|n-1, \theta) + \theta \cdot \text{Bin}(x-1|n-1, \theta) \Big] \\[6pt]
&= (1-\theta) \times \sum_{x=0}^{n-1} f_x \cdot \text{Bin}(x|n-1, \theta) + \theta \times \sum_{x=0}^{n-1} f_{x+1} \cdot \text{Bin}(x|n-1, \theta) \\[6pt]
&= (1-\theta) \cdot B(\triangleleft \ \mathbf{f}, \theta) + \theta \cdot B(\triangleright \ \mathbf{f}, \theta), \\[6pt]
\end{aligned} \end{equation}$$
donde $\triangleleft \ \mathbf{f} = (f_0,...,f_{n-1})$ e $\triangleright \ \mathbf{f} = (f_1,...,f_{n})$. Esto nos da una ecuación recursiva de la curva Bézier, donde la recursividad se descompone la función en una suma ponderada de la curva Bézier para los más pequeños de control de vectores. En cada paso de la iteración, la longitud del vector de control se reduce en un elemento.
De Casteljau del algoritmo: Este algoritmo tiene la ventaja de la anterior ecuación recursiva de la curva Bézier. Para explicar el algoritmo, se definen los operadores de $\triangleleft$ e $\triangleright$ a truncar el argumento del vector por un elemento de la derecha y la izquierda respectivamente. Ahora, teniendo un valor fijo de $\theta$ podemos definir los valores:
$$B_{k,x} \equiv B(\triangleleft^{n-k-x} \triangleright^x \mathbf{f},\theta),$$
y podemos organizar estos valores en un $(n+1) \times (n+1)$ matriz $\mathbf{B}$ como sigue:
$$\mathbf{B} \equiv \begin{bmatrix}
B_{0,0} & B_{0,1} & B_{0,2} & \cdots & B_{0,n-2} & B_{0,n-1} & B_{0,n} \\
B_{1,0} & B_{1,1} & B_{1,2} & \cdots & B_{1,n-2} & B_{1,n-1} & 0 \\
B_{2,0} & B_{2,1} & B_{2,2} & \cdots & B_{2,n-2} & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
B_{n-2,0} & B_{n-2,1} & B_{n-2,2} & \cdots & 0 & 0 & 0 \\
B_{n-1,0} & B_{n-1,1} & 0 & \cdots & 0 & 0 & 0 \\
B_{n,0} & 0 & 0 & \cdots & 0 & 0 & 0 \\
\end{bmatrix}.$$
La parte inferior del elemento de esta matriz es $B_{n,0} = B(\triangleleft^0 \triangleright^0 \mathbf{f},\theta) = B_\mathbf{f}(\theta)$, que es la función de salida que queremos calcular, y la fila superior se compone de los valores de $B_{0,x} = B(\triangleleft^{n-x} \triangleright^x \mathbf{f},\theta) = f_x$, que son las iniciales de "puntos de control" de la función.
De Casteljau del algoritmo comienza con el vector de control en la parte superior de la línea de esta matriz, y trabaja hacia abajo a través de la anterior matriz de valores mediante la anterior recursiva caracterización de la curva Bézier. En escalares de la forma, el recurrente ecuaciones para el algoritmo son:
$$\begin{equation} \begin{aligned}
B_{0,x} &\equiv f_x
\quad \quad \quad \quad \quad \quad \quad
\quad \quad \quad \quad \quad \quad \quad \text{ }
\text{ for } x = 0,...,n, \\[6pt]
B_{k,x} &\equiv (1-\theta) \cdot B_{k-1,x} + \theta \cdot B_{k-1,x+1}
\quad \quad \quad \text{for } x = 0,...,n-k.
\end{aligned} \end{equation}$$
Como se puede observar en la matriz, este algoritmo consiste en el cálculo de $n(n+1)/2$ valores del conjunto inicial de puntos de control, por lo que tiene la complejidad de la $\mathcal{O}(n^2)$. Cada recursiva de cálculo es un simple promedio ponderado de los elementos de arriba y de arriba-derecha, de manera que la carga computacional de cada paso implica dos multiplicaciones y una adición (y por lo tanto es relativamente pequeño).
El algoritmo puede ser implementado en el estándar de la probabilidad de escala, siempre y cuando el ambiente computacional es suficiente para evitar subdesbordamiento de problemas para las pequeñas probabilidades. Alternativamente, el cálculo se puede hacer en log-probabilidad de escala para evitar subdesbordamiento de problemas (véase aquí para una discusión acerca de la adición de pequeñas probabilidades en el registro de la escala).
