Linear fit to data samples from a logistic function Assume that I uniformly sample data points from an interval [-$a$; $a$], then pass these data points to a logistic function of the form
$f(x) = \frac{1}{1+e^{-\beta\cdot x}}$,
and fit a linear line to these data samples (x, f(x)). Can the expected slope of this linear fit be expressed as a function of $\beta$ and $a$?
(Bonus: tangens hyperbolicus function $f(x)=\tanh(\beta x)$ instead of the logistic function)
 A: You have data $X_1,\dots,X_n\stackrel{\text{iid}}\sim\text{Unif}(-a,a)$ and obtain $Y_i = \sigma(X_i)$ where $\sigma(z) = \frac{1}{1+e^{-bx}}$.
The population slope of a simple linear regression is$\newcommand{\E}{\operatorname{E}}$
$$
\beta := \frac{\text{Cov}(X_i, Y_i)}{\text{Var}(X_i)} = \frac{\E(X\sigma(X)) - (\E X) (\E \sigma(X))}{\text{Var}(X)}.
$$
$X$ is symmetric about $0$ and the variance is just that of a uniform RV which is known, so all we really need to compute is $\E(X\sigma(X))$. This is
$$
\E(X\sigma(X)) = \frac 1{2a} \int_{-a}^a \frac{x}{1+e^{-bx}}\,\text dx.
$$
We can note that
$$
\int \frac{1}{e^{-bx} + 1}\,\text dx = \int \frac{e^{bx}}{1 + e^{bx}}\,\text dx \\
= \frac 1b \int \frac 1u \,\text du = \frac 1b \log (e^{bx} + 1)
$$
so we can integrate by parts to get
$$
\E(X\sigma(X)) = \frac{1}{2ab}x\log(e^{bx}+1)\bigg\vert_{-a}^a - \frac 1{2ab}\int_{-a}^a \log(e^{bx} + 1)\,\text dx.
$$
With the first term (aside from some scaling constants) we end up with
$$
\log(e^{ab}+1) + \log(e^{-ab}+1) = \log\left[e^{ab}(e^{-ab}+1)\right] + \log(e^{-ab}+1) \\
= ab + 2\log(e^{-ab}+1)
$$
so
$$
\frac{1}{2ab}x\log(e^{bx}+1)\bigg\vert_{-a}^a = \frac a2 + \frac 1b \log(e^{-ab}+1).
$$
For the other term, we can let $u = -e^{bx}$ so $\frac 1{bu}\,\text du = \text dx$ which means$\newcommand{\Li}{\operatorname{Li}_2}$
$$
\int_{-a}^a \log(e^{bx} + 1)\,\text dx = -\frac 1b \int_{-e^{ab}}^{e^{ab}} \frac{\log(1-u)}u\,\text du \\
= \frac 1b\left[\Li(e^{ab}) - \Li(-e^{ab})\right]
$$
where $\Li$ is the dilogarithm function. All together this means
$$
\E(X \sigma(X)) = \frac a2 + \frac 1b \log(e^{-ab}+1) + \frac{\Li(-e^{ab}) - \Li(-e^{-ab})}{2ab^2}.
$$
I have something of the form $\Li(z) - \Li(1/z)$. Using the fact that
$$
\Li(z) + \Li(1/z) = -\frac{\pi^2}6 - \frac 12 \log^2(-z)
$$
I can rewrite this with just a single $\Li$ in it:
$$
\E(X \sigma(X)) = \frac a2 + \frac 1b \log(e^{-ab}+1) + \frac{-\frac{\pi^2}6 - \frac 12 a^2b^2 - 2 \Li(-e^{-ab})}{2ab^2} \\
= \frac a4 + \frac 1b \log(e^{-ab}+1) - \frac{\pi^2 + 12 \Li(-e^{-ab})}{12ab^2}.
$$
This means
$$
\beta(a,b) = \frac{\frac a4 + \frac 1b \log(e^{-ab}+1) - \frac{\pi^2 + 12 \Li(-e^{-ab})}{12ab^2}}{a^2/3} \\
= \frac{3}{4a} + \frac{3}{a^2b}\log(e^{-ab}+1) - \frac{\pi^2 + 12 \Li(-e^{-ab})}{4a^3b^2}.
$$
$\Li$ is a special function and can't be expressed in terms of elementary functions in general, and since this could be solved for $\Li$ that means there's no general elementary expression for $\beta$. But we can evaluate it in terms of elementary functions for particular values of $a$ and $b$. One such value is
$$
\Li\left(-\varphi\right) = -\frac{\pi^2}{10} - \log^2\varphi
$$
so if $ab = -\log(\varphi)$ then we'll be able to evaluate $\beta$ in terms of elementary functions. As an example of this, we could take
$$
a = \sqrt 2 \\
b = -\frac 1{\sqrt 2}\log(\varphi).
$$
Although this is still a really messy expression even if it's only in terms of elementary functions.
This was a lot of work just to express the integral $\int_{-a}^a \frac{x}{1+e^{-bx}}\,\text dx$ in terms of a different integral, but I think the value is that it shows that there isn't a closed form for $\beta(a,b)$ that we're missing, and it relates it to a well-studied special function that has high quality implementations available.
Here's a simulation to check.
set.seed(111)
nsim <- 1e6
a <- 2.34
b <- 1.2
x <- runif(nsim, -a, a)
y <- 1 / (1 + exp(-b * x))

curve(plogis(b*x), -a, a, 500, col=4)
points(y[1:100]~x[1:100], cex=.5)

# I'm just integrating to avoid needing to load other libraries
Li2 <- function(z) -integrate(function(u) log(1-u)/u, 0, z)$value
slopefunc <- function(a,b) {
  3/(4*a) + 3 / (a^2 * b) * log(exp(-a*b) + 1) - (pi^2 + 12*Li2(-exp(-a*b))) / (4*a^3*b^2)
}

lm(y~x)$coef[2]
slopefunc(a,b)
cov(x, y) / var(x)

