Line of best fit (Linear regression) over vertical line I want to get a line of the best fit which is a line that passes as close as possible to a set of points defined by coordinates point_i = (X_i, Y_i).
When I apply linear regression, I have a special case where the line do not correspond to what I want. It is the case where the data-points are aligned over a vertical line (i.e. Y values do not depend on X values) like the following example for instance (or see the points on image bellow), e.g:
Y = [1,2,3,4,5,8,10,11,15,13,21,22,23,24,25,28,30,31,35,33]
X = [1,1,1,1,1,1,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ];

In this case I expect to get a line of best fit which is vertical and passes as close to a set of points, but when I get from regression is a horizontal line ! However, what I'm interested in, is not the error over the Y axis, but I'm interested in the error to minimize which indicate how far (according to distance) are the points from the line of best fit that we get (e.g. the sum of shortest distance from each point to the line).
Is there any why to get such a line and to get the corresponding score which indicates how far are the points this line ?
I think "linear regression" is not the same thing as "line of best fit", and what I want to get is "line of best fit".

EDIT:
I want to estimate line of best fit for another purpose, not to evaluate how much X and Y are correlated (it is more for a pattern recognition problem, where my data-points are some pixels defined by their coordinates x,y on the image). So I want to fit the line even if X values are not correlated to Y values. And if I just always swap X and Y to be able to fit the line, the next time if I have points lying on a horizontal line (Y = [3,3,3,3], X=[3, 6, 23, 30] for example) the if they are swapped I'll get again the same problem. I want to always fit the line whatever is the correlation between W and Y (even if there is no correlation, like in the example I provided).
So to be more clear: Whatever is the correlation between X and Y, my objective is not to predict some y values based on new values of x, my objective is rather just to get a line which passes as close as possible to all points that we have (even if Y are independent of X ...).
 A: One way to do this would be to check the point "verticality", regress a line through (X,Y) inverted points, and then invert the regressed line back.


*

*Check point "verticality" by comparing the standard deviations of X and Y values.

stDevX = StandardDeviation(xValues) 
stDevY = StandardDeviation(yValues) 

verticality = stDevY/stDevX #Handle div/0 cases




*If vertical, regress a line through inverted Y and X points.

isVertical = verticality > 1.5 #adjust based on situation (could be 2.0 or more)

if(isVertical)
   invertedLine = RegressLine(yValues,xValues) #Note reversed x,y order




*Obtain non-inverted slope and intercept from inverted line

slope = 1/invertedLine.Slope
intercept = -invertedLine.Intercept / invertedLine.Slope


Notes: 


*

*This will fail if the points are perfectly vertical (ie. stDevX is
=0). If this is the case, the line equation is just x=Average(xValues). For nearly-vertical lines, the slope will be
large, but finite.

*The stDevY/stDevX ratio will be approximately equal to the regressed line |slope|. When the points become vertical, the regressed line slope will start to deviate from the ratio. In my case, whenever the ratio exceeded 1.5, the regressed slope became inaccurate, and flipping the axes improved accuracy. You should experiment with your data to find the best verticality threshold.

A: In the case you want to make a linear regression where Y values do not depend on X values, the method is to look for the solution of equation ny + p = x instead of the usual mx + p = y

The principle is the same, using the least square method. It solves the equation ax = b by computing a vector x that minimizes the Euclidean 2-norm || b - a x ||^2
I've made a script where I compare the two methods (ny + p = x versus mx + p = y) using 3 lines, a diagonal, a horizontal line and a vertical one. I'm displaying the Sum of the residuals (Sum of Res.) in the legend.
You could see that for the diagonal line, the results are similar, but for the vertical and horizontal one, 1 method is giving good results while the other is giving is not.
Method 1: mx + p = y
The green line equation is not well fitting the values. The sum of residuals is very high compared to the good regression fits like diagonal or horizontal lines.

Method 2: ny + p = x
The red line equation is not well fitting the values. The sum of residuals is very high compared to the good regression fits like diagonal or vertical lines.

A: One purpose of a linear model is to see how a variable varies as a function of variation in another variable - it's basically just a re-parameterized covariance. In your example, you have a variable with zero variation, thus it would be nonsensical to use linear regression. If you still want to fit the line, you could just swap X and Y. The resultant regression coefficient would of course still be zero (indicating no relationship, or slope), and further confirm that knowledge of X (using your original notation) tells us nothing about Y. Given this is true, the question you should answer is, "Why do I want to estimate a line of best fit?" 
A: The reason why this does not work is because your system (or data, rather) does not fit the linear model.
Your data vector had the problem that each element is exactly the same. In other words if you'd have multiple data vectors you'd have a linear dependend system.
Or put differently yet again - the issue here is perfect multicollinearity.  
This is a situation where the linear model is actually not able to give you any information. It might be that you didn't see this point when studying this because it usually comes in a form of "Matrix X is of rank N" or something like that, but it still important.
Think about it:
If your regressand is always 1 as indicated in your data set, it is obviously impossible to put any line through it as the variation of y can not be related to a constant factor. The best approximation, linear wise, is the mean of $y$ which is what you seem to get. This is still the closest you can functionally get to $y$.  
But once again just so that you can understand this: If your best functional approximation is a horizontal line, a mean, then the actual data points have to be created by an error term. Intuitively this makes sense: If your x is always the same, how would you ever attribute influence of it on a changing term?  
In the end you are looking - in your graph at least - for a measure based on only $x$. In your graph this would be the arithmetic mean of $x$. You can see that even in your graph the $y$ have no bearing on the positioning of your $y$.
A possible solution to your question is actually to flip $x$ and $y$.
What you are looking for then are the residuals of that regression as a measure of "closeness" to the line. This is assuming your $x$ is not actually equal to 1 everywhere.
Use this:  $$x=\alpha+y+\epsilon$$ 
If you are fixed on getting a vertical straight line (not a function) use the model without regressand instead of which will give you the mean of $x$ instead as " horizontal line": $$x=\alpha+\epsilon$$
If you flip that graph, you have your vertical line - or your fitted "almost horizontal line" if you used the first model.
Next, calculate your residuals. These happen to be the differences between your $y$ values and your predicted $y$ values - the fitted or vertical line depending on the model.
This is essentially what you want I guess, though honestly I don't know why.
You can now calculate the sum or residuals, their squares and for that the standard errors and such...
Edit: I am going to assume your edit was targeted at my post?
If you happen to have this problem in both sides you'll have to catch those special cases anyway. Regression technically can not be done on this fringe case as I posted above. You can see this even because you had to modify your graph to display variation in $x$ even though your data says the $x$ is always 1.
In practice, why would you get cases where entries are literally the same? The linear model stays accurate even for cases with a lot of multicollinearity (you don't have to flip n that case).  
My final tip would be to add random jitter to your $x$ or $y$ so you don't get that exact problem. But if you are really looking at cases where you can not catch and remove cases of perfect multicollinerity, you are probably asking the wrong question or model tbh
