I have come upon with the next expression:
\begin{equation} P_r=\prod_{k=0}^r \left(1+\frac{1}{2}\left(\frac{\frac{1}{2}+k+1}{\frac{1}{2}+k} -\frac{\frac{1}{2}+k}{\frac{1}{2}+k+1}\right)\right) \end{equation}
And simulations have shown that its value can be approximated quite accurately to a linear function in the form: \begin{equation} P_r \approx O+M \cdot \left(\frac{1}{2}+r \right) \end{equation}
Using such knowledge, I have tried to obtain values for $O$ and $M$ using the first two values of $P_r$, that is, for $r=0$ and $r=1$. This approach is only reasonably accurate for small values of $r$, but tends to diverge from the actual value for higher values of $r$.
Could anyone help finding out the theoretical derivation to obtain a good linear approximation, and therefore the most accurate values for $O$ and $M$?
Some Matlab code:
x = [5];
x0 = 5;
x_r = zeros(200,1);
x_r(1) = x0/2;
P_r = zeros(201,1);
P_r(1) = 1;
x_r(2) = x(randi(length(x)));
zeta_r(2) = sum(x_r(1:2));
P_r(1) = (1+1/2*(sum(x_r(1:2))/sum(x_r(1:1))-sum(x_r(1:1))/sum(x_r(1:2))));
for r=3:201
x_r(r) = x(randi(length(x)));
zeta_r(r) = sum(x_r(1:r));
P_r(r-1) = P_r(r-2)*(1+1/2*(sum(x_r(1:r))/sum(x_r(1:r-1))-sum(x_r(1:r-1))/sum(x_r(1:r))));
end
% P_r = P_r(2:201);
hold on;
plot(cumsum(x_r(1:200)),P_r(1:200));
%Using two first points
M = (P_r(2)-P_r(1))/x_r(2);
O = P_r(1)-M*x0/2;
plot(cumsum(x_r),O+M*cumsum(x_r),'r');
%Using first and third points
M = (P_r(3)-P_r(1))/(x_r(2)+x_r(3));
O = P_r(1)-M*x0/2;
plot(cumsum(x_r),O+M*cumsum(x_r),'black');
%Using first and some other point
M = (P_r(10)-P_r(1))/(sum(x_r(2:10)));
O = P_r(1)-M*x0/2;
plot(cumsum(x_r),O+M*cumsum(x_r),'green');
h = fit(cumsum(x_r(1:200)),P_r(1:200),'poly1');
plot(h,cumsum(x_r(1:200)),P_r(1:200));
Hint. One may first observe that $$ 1+\frac{1}{2}\left(\frac{1/2+k+1}{1/2+k}-\frac{1/2+k}{1/2+k+1}\right)=\frac{\left(k+\frac{3+\sqrt{2}}2\right) \left(k+\frac{3-\sqrt{2}}2\right)}{\left(k+\frac12\right)\left(k+\frac32\right)} \tag1 $$ then, using the $\Gamma$ function, we get $$ \prod_{k=0}^{r}(k+a)=\frac{\Gamma(r+a+1)}{\Gamma(a)} \tag2 $$ applying it to $(1)$ gives the following closed form of the initial product $$ \begin{align} &\prod_{k=0}^{r}\left(1+\frac{1}{2}\left(\frac{1/2+k+1}{1/2+k}-\frac{1/2+k}{1/2+k+1}\right)\right) \\\\&=\frac{\Gamma\left(r+\frac{5+\sqrt{2}}2\right) \Gamma\left(r+\frac{5-\sqrt{2}}2\right)\Gamma\left(\frac12\right)\Gamma\left(\frac32\right)}{\Gamma\left(\frac{3+\sqrt{2}}2\right)\Gamma\left(\frac{3-\sqrt{2}}2\right)\Gamma\left(r+\frac32\right)\Gamma\left(r+\frac52\right)}, \quad r\geq0.\tag3 \end{align} $$
or, with some simplifications,
The function $r \mapsto f(r)$ on the right hand side of $(4)$ is $C^{\infty}$ over $r \in [0,\infty)$. It then admits a linear approximation on the form $ f(r)=f(r_0)+f'(r_0)(r-r_0)$ near each $r_0 \in [0,\infty)$.
For example, as $r$ is small, one finds that
and, using the asymptotic expansion of the $\Gamma $ function, as the argument $r$ is great, one gets