def sample_mean(xs):
    return sum(xs) / len(xs)

def sample_std(xs):
    m = sample_mean(xs)
    return (sum([(i-m)**2 for i in xs]) / (len(xs)-1))**0.5

def sample_cov(xs, ys):
    x_mean = sample_mean(xs)
    y_mean = sample_mean(ys)
    o = 0.0
    for index in range(len(xs)):
        o += (xs[index] - x_mean)*(ys[index] - y_mean)
    return o / (len(xs) - 1)

def pearson_r(xs, ys):
    return sample_cov(xs, ys) / (sample_std(xs)*sample_std(ys))

def linear_regression(xs, ys):
    r = pearson_r(xs, ys)
    m = r * sample_std(ys) / sample_std(xs)
    b = sample_mean(ys) - m * sample_mean(xs)
    return m, b, r