Fit a linear function to multiple measurements I have the data of a measurement of the same value that was repeated multiple times to decrease random noise. 
There are multiple values per input-value (time), an example could look like this:
Time|Measurement 1|Measurement 2|Measurement 3|Measurement 4 ...
x_t |y_t1         |y_t2         |y_t3         |y_t4
0   |value        |value        |value        |value ...
1   |value        |value        |value        |value ...
2   |value        |value        |value        |value ...
3
...
...


How can I find the linear function $y_t=\beta x_t +c$ that fits "best" to the data? Also, how can I then calculate the standard deviation/error on the slope and intercept?
One approach I thought of, is to just take the average of all the measurements and then perform linear regression on the average measurement but that's probably not the best method.
 A: I used the analytical method of linear regression and regressed over the data as if it was just one dataset (as proposed by James Phillips' comment on the question). I don't know if I should even post my code here since this is not stackoverflow, but it might help someone someday:
import numpy as np

def mult_linreg(x, y):
    if len(x) < 3:
        raise Exception("More than two datapoints required")
    sum_xy = 0
    n = 0
    sum_y = 0
    sum_x = sum(x) * len(y)
    sum_x_sq = sum(x**2) * len(y)

    for y_ in y:
        if len(y_) is not len(x):
            raise Exception("x and y dimension mismatch")
        sum_xy += np.dot(x, y_)
        n += len(y_)
        sum_y += sum(y_)
    nxsq_xsq = (n*sum_x_sq - sum_x**2)
    a = (n * sum_xy - sum_x * sum_y)/nxsq_xsq
    b = (sum_x_sq * sum_y - sum_x * sum_xy)/nxsq_xsq
    sum_y_ax_b_sq = 0
    for i in range(0, len(x)):
        xi = x[i]
        axi_and_b = a*xi + b
        for y_ in y:
            yi_ = y_[i]
            sum_y_ax_b_sq += (yi_ - axi_and_b)**2

    S = np.sqrt(sum_y_ax_b_sq/(n-2))
    err_a = S * np.sqrt(n/nxsq_xsq)
    err_b = S * np.sqrt(sum_x_sq/nxsq_xsq)
    return a,b, err_a, err_b

#example
print(mult_linreg(np.array([0,6,10]),np.array([[-1,2,4],[-1.1,2.1,3.9]])))

The result is indeed very similar to what you get when regressing the average of all datasets, but that might depend on the input
