Calculation of aX+bY standard deviation and covariance of X and Y

[Authenian]

def randcal(X,Y,a,b):
    import numpy as np
    from math import sqrt
    X = np.array(X)
    Y = np.array(Y)
    EX = np.mean(X)
    EY = np.mean(Y)
    Z = X*Y
    EZ = np.mean(Z)
    COV = EZ - EX*EY
    if COV > 1e-3 or COV < -1e-3:
        print("X and Y are dependent")
    X2 = X**2
    Y2 = Y**2
    DX = np.mean(X2) - EX**2
    DY = np.mean(Y2) - EY**2
    sigmaX = sqrt(DX)
    sigmaY = sqrt(DY)
    print(f"sigma X = {sigmaX}")
    print(f"sigma Y = {sigmaY}")
    DXY = a**2*DX + b**2*DY + 2*a*b*COV
    sXY = sqrt(DXY)
    print(f"sigma aX+bY = {sXY}")
    return COV

[Authenian]

Use this for example: randcal([3, 6, 2, 8, 4],[2, 1, 9, 6, 0],2,2)