linear regression, symmetry of model does not lead to symmetry of coefficients Experiment: You are given a large population of real numbers. For simplicity take the whole numbers from -n to n. Take two independent random samples x and y of size k and sort them (each one individually). Compute the linear regression $y=a+bx$. 
Question: what are the expected values for a and b? What does the distribution look like?
As the setting is symmetric in x and y I expected $a=0$ and $b=1$ but numerics show that this is not the case.
I first posted this question in the math section:
https://math.stackexchange.com/questions/2476502/linear-regression-symmetry-of-model-does-not-lead-to-symmetry-of-coefficients?
but did not get any useful answers. This question here
symmetry of linear regression
is also closely related.
Here is some R code that runs this experiment. 
n_pop  <- 100

n_sample <- 20

draws <- 1000

M <- matrix(c(1:(2*draws)), nrow=draws)

for (j in 1:draws) {

  x <- sort(sample(-n_pop:n_pop,n_sample,replace=FALSE))

  y <- sort(sample(-n_pop:n_pop,n_sample,replace=FALSE))

  xy_model <- lm(y~x)

  M[j,] <- coef(xy_model)

}

mean(M[,1])

mean(M[,2])

sd(M[,1])

sd(M[,2])

hist(M[,1],breaks=20)

hist(M[,2], breaks=20)

t.test(M[,2], mu=1, conf.level=0.99)

t.test(M[,1], mu=0, conf.level=0.99)

For this population and sample size I get around $a=1$ and $b=0.98$, $a=0$ and $b=1$ can be excluded with high confidence.
I suspect this is related to the fact that OLS minizes distance only on the y-axis but I don't understand why this changes the coefficients. Numerically, if the sample size gets smaller, the cofficients move farther way from 0 and 1.
If one uses a Deming regression which minimizes euclidean distance to the line, numerically $a=0$ and $b=1$ looks plausible. 
edit: the code above gives a distribution for $a$ that is symmetric around 0, I had an error in my R code. The question for $b$ remains and the answer below gives an excellent explanation why $b$ is not equal to $1$.
 A: 
1
The regression $y = a + bx + \epsilon$ is not symmetric for replacing y and x, as you already suggested via the link (symmetry of linear regression).
However, it is not trivial whether the expectation values, $\hat{a}$ and $\hat{b}$, are going to be different from a symmetric $E(\hat{a}_{xy})=E(\hat{a}_{yx})=0$ and $E(\hat{b}_{xy})=1/E(\hat{b}_{yx})=1$. Possibly the difference in the regression $y \sim x$ vs $x \sim y$ cancel over all possibilities (we will see this is true for $a$). For this we need to "see" all the possibilities and how they average/cancel.

2
We will show that the intercept is symmetric
$$E(\hat{a})=0$$
however the slope is not, and will be smaller than 1
$$E(\hat{b})<1$$
We show this by comparing the non-symmetric regression with a symmetric regression. 

3
As symmetric regression we use the line in between the two regression lines lines defined by the y~x and x~y.
In algebraic terms:
$$\hat{b}_{sym} = \frac{\hat{b}_{xy} + \hat{b}_{yx}^{-1}}{2}$$
and 
$$\hat{a}_{sym} = \frac{\hat{a}_{xy} - \hat{a}_{yx}\hat{b}_{yx}}{2}$$

4
Note the following about the symmetric slope. 


*

*The slope will be symmetric around the angle of 45 degrees 
ie, a probability for an angle with $45+\alpha$ is equal to a
probabilty for an angle with $45-\alpha$, also the expectation value
for its angle coefficient will be 1, since:


*

*$\hat{b}_{sym}(x,y)=\hat{b}_{sym}(y,x)$, thus
$E(\hat{b}_{sym}(x,y))=E(\hat{b}_{sym}(y,x))$ 
(note that we do not have the same for the slope of the asymmetric case $\hat{b}_{xy}(z_1,z_2)=\hat{b}_{yx}(z_2,z_1)$ ) 

*yet we must also have  the symmetry $E(\hat{b}_{sym}(x,y))=1/E(\hat{b}_{sym}(y,x))$, the line that is described by the expectation of the coefficient, if we change the labels xy than we expect a mirror image which has a different b coefficient (the inverse). 


*The slope $\hat{b}_{xy}$ will be smaller due to regression to the mean.


Therefore... the symetric slope will be $E(\hat{b}_{sym})=1$, however then since always $\hat{b}_{xy} \leq \hat{b}_{sym}$ we will have $E(\hat{b}_{xy}) \leq E(\hat{b}_{sym})=1$. And the equality is only true when there are no error terms.
See the following image to recall the regression lines $xy$, $yx$ (red curves) and the symmetric one in between (green curve), this image is from the example code below: 

If we plot the angle of the symmetric line (green) and the regression line for the y~x model (the red one that has smaller angel with the x-axis) then we see that the angle of the symmetric line distributes evenly around 45 degrees, however the angle of the non-symmetric line  $y~x$, is always lower: 


5
What is left is to show that $E(\hat{a})=0$. I do not have a direct proof for this. However, it should be sufficient to notice that the regression to the mean (which is zero) does not favor either larger or smaller $a$. Note that there is a point-symmetry, $(x,y) \rightarrow (-x,-y)$, for which the angle $b$ is invariant, but not the intercept $\hat{a}$, which means that it, $a$, should be zero.
You obtain values of $\hat{a}$ different form zero, but you should take into account it's sample variance. (The sample error for $\hat{a}$ is different from zero)

6
Below is adjustment to your code to examine the above described effects
n_pop  <- 100

n_sample <- 20

draws <- 1000

M <- matrix(c(1:(2*draws)), nrow=draws)
K <- matrix(c(1:(2*draws)), nrow=draws)
ma <- matrix(c(1:(2*draws)), nrow=draws)
meanx <- matrix(c(1:(1*draws)), nrow=draws)

for (j in 1:draws) {

  x <- sort(sample(-n_pop:n_pop,n_sample,replace=FALSE))

  y <- sort(sample(-n_pop:n_pop,n_sample,replace=FALSE))

  xy_model <- lm(y~x)
  yx_model <- lm(x~y)

  M[j,] <- coef(xy_model)      
  K[j,] <- coef(yx_model)     

  # xy and yx model coefficients last loop
  b1 <- coef(xy_model)[2]
  b2 <- 1/coef(yx_model)[2]
  a1 <- coef(xy_model)[1]
  a2 <- -coef(yx_model)[1]*coef(yx_model)[2]

  ma[j,1] <- 0.5*sum(a1,a2)
  ma[j,2] <- 0.5*sum(b1,b2)

  meanx[j,1] <- mean(x)
}

# means
mean(M[,1])
mean(M[,2])
mean(ma[,1])
mean(ma[,2])



# plotting last for-loop x-y as example of regresion fits
plot(x,y)

xp <- -n_pop:n_pop
lines(xp,a1+b1*xp,col=2)
lines(xp,a2+b2*xp,col=2)
lines(xp,0.5*a2+0.5*a1+0.5*b2*xp+0.5*b1*xp,col=3)


# plotting coefficients symmetric versus xy
plot(180/pi*atan(ma[,2]),180/pi*atan(M[,2]), xlab="symmetric fit angle", ylab="xy-fit angle")
lines(c(0,90),c(0,90),col=2)

# plotting mean of x vs beta
plot(meanx,ma[,2]-M[,2], xlab="mu_x",ylab="beta_sym - beta_xy",pch=21,
     col=adjustcolor("black",alpha.f=0.1),
     bg=adjustcolor("black",alpha.f=0.1),log="y")

