Regression when output is in a specific interval Say my output is continuous on the interval $[-1, 1]$. Is there a standard technique to deal with this? Do I just use logistic function, same as in logistic regression, but for actual regression rather than classification? What are other options, what is the standard approach?
 A: The simple linear regression theory is more developed for normal variables than for other distributions. When we have to deal with a problem like yours, we can use change of variables. In your case, I will use a change like:
$$ z = \frac {2y} {1-y^2} = \frac {1} {1-y} - \frac {1} {1+y} $$
This function is increasing: if $y$ is greater, $z$ is greater. When $y$ is near $-1$, $z$ is near $-\infty$; when $y$ is near $+1$, $z$ is near $+\infty$. With these trick, you can calculate the linear relationship between the independent variable $x$ and the the dependent variable $z$.
A: The appropriate technique depends on your goal.
If you are building a model for inference, you should focus on the properties of the distribution of your target conditional on covariates, $p(y|x)$. 
For example, the value $0.5(y+1)$ may be distributed as $Beta(\alpha(x), \beta(x))$. In this case, you may perform maximum likelihood estimation of parameters of the functions $\alpha(x)$ and $\beta(x)$, and find the best form for them (e.g. linear or log-linear). Google "beta regression" for more details.
Instead of $Beta$, you can fit a GLM with any link function you want (indeed, logit link is commonly used). You can also map $y$ into $(-\infty, \infty)$ with any function you want, and use unconstrained regression. The last approach, however, can fail if exact $\pm 1$s are present in your data.
Another trick is to transform your regression into weighted classification. From each training observation $(x, y)$ you can generate two observations $(x, 1)$ and $(x, 0)$ with corresponding weights $\frac{1+y}{2}$ and $\frac{1-y}{2}$, fit a probabilistic classifier (e.g. logistic or probit regression), and then transform predicted probability of $1$ back to $y$.
If you are building a model for prediction, the probabilistic properties may be ignored, you just focus on predicting $y$ as close as possibly, whatever it means. In this case, you may fit any function $y=f(x)$, and just truncate outside $[-1, 1]$. This approach allows you to try lots of different regression algorithms without bothering much about the boundaries on $y$. 
Moreover, several machine learning models (e.g. decision trees and their ensembles random forests, k-nearest-neighbor, or any other method which prediction a is weighted average of training samples) are by design unable to predict higher than the highest training value, or lower than the lowest. If you use them, you may never worry about the interval of $y$.
What approach is standard, depends on the domain and on your goal. But fitting a logistic function to continuous data seems to be OK:


*

*it always predicts in $(-1, 1)$

*it works even with exact $\pm 1$

*generalized linear form gives you a basis for inference and feature selection

*it had decent prediction accuracy in the most cases I seen.


Now it's time for an implementation. There is an example of R code that evaluates such a model. 
set.seed(1)
data = data.frame(x=1:100)
data$y = 1 / (1 + exp(5-0.1*(data$x) + rnorm(100)))

model = glm(y~x, family = 'binomial', data=data)
summary(model)
plot(x, y)
lines(x, predict(model, data, type = 'response'))

It outputs the following table of estimated coefficients (close to the "true" coefficients I used)
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.48814    0.88243  -5.086 3.65e-07 ***
x            0.08713    0.01615   5.394 6.89e-08 ***

and a picture with the training data and the fitted function

Unfortunately, Python's sklearn does not allow logistic regression to run in regression mode, but it is possible with statsmodels - it has a Logit class that allows continuous targets. The interface and output are pretty similar to those in R:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
np.random.seed(1)
df = pd.DataFrame({'x': range(100)})
df['y'] = 1 / (1 + np.exp(5-0.1*(df.x) + np.random.normal(size=100)))
model = smf.logit('y~x', data=df).fit()
print(model.params)
plt.scatter(df.x, df.y)
plt.plot(df.x, model.predict(df), color='k')
plt.show()

One more issue worth considering is evaluation metric for your model. Along with standard RMSE and MAE, in such problem rank-based metrics, such as Spearman correlation, may be useful. If you do weighted classification instead of regression, you can also calculate weithted classification metrics, like ROC AUC. 
The rationale for such metrics is that in the end you may want not to predict $y$ as accurately as possible, but separate low $y$ from high $y$ as accurately as possible, but you don't know the threshold in advance, or it is variable. Rank-based metrics reflect this process better than difference-based metrics.
