【測繪程式設計試題集】 試題03 利用線性迴歸模型進行衛星軌道的預報

資料

300,21182.88,-7044.56,14639.48
600,21707.87,-6930.28,13906.68
900,22207.04,-6828.65,13147.66
1200,22679.16,-6738.66,12363.84
1500,23123.06,-6659.23,11556.71
1800,23537.69,-6589.21,10727.78
2100,23922.07,-6527.40,9878.61
2400,24275.33,-6472.54,9010.81
2700,24596.67,-6423.32,8126.00
3000,24885.42,-6378.40,7225.86
3300,25141.01,-6336.41,6312.08
3600,25362.96,-6295.93,5386.38
3900,25550.92,-6255.54,4450.51

問題

解

def read_m(path):
    # 　所有資料
    m = []
    # x
    xlist = []
    # y
    ylist = []
    # z
    zlist = []
    # time
    time_list = []

    with open(path, 'r') as f:
        for i in f.readlines():
            aa = i.replace('\n', '').split(",")
            bb = [eval(a) for a in aa]
            m.append(bb)
            time_list.append(bb[0])
            xlist.append(bb[1])
            ylist.append(bb[2])
            zlist.append(bb[3])
    return {
        "alldata": m,
        "time": time_list,
        "x": xlist,
        "y": ylist,
        "z": zlist,
    }


XXX = None
YYY = None


def xpj():
    """
    X平均值
    :return:
    """
    sum = 0

    for i in range(XXX.__len__()):
        sum += XXX[i]

    return sum / XXX.__len__()


def ypj():
    """
    Y 平均值
    :return:
    """
    sum = 0

    for i in range(YYY.__len__()):
        sum += YYY[i]

    return sum / YYY.__len__()


def sse():
    """
    迴歸方程
    :return:
    """
    sum = 0
    xa = xpj()
    ya = ypj()

    for i in range(XXX.__len__()):
        sum += (XXX[i] - xa) * (YYY[i] - ya)

    return sum


def ssx():
    """
    迴歸方程
    :return:
    """
    sum = 0
    xa = xpj()
    for i in range(XXX.__len__()):
        sum += (XXX[i] - xa) * (XXX[i] - xa)
    return sum


def getbeta1():
    """
    bate1
    :return:
    """
    bbeta = sse() / ssx()
    return bbeta


def getbeta0():
    """
    beta0
    :return:
    """
    return ypj() - getbeta1() * xpj()


def huiguixishu(x, y):
    """
    迴歸係數
    :param x:
    :param y:
    :return:
    """
    global XXX
    global YYY
    XXX = x
    YYY = y

    beta1 = getbeta1()
    beta0 = getbeta0()
    return [beta0, beta1]


def predic(x, beta0, beta1):
    """
    估計
    :param x:
    :param beta0:
    :param beta1:
    :return:
    """
    a = beta0 + beta1 * x
    return a


if __name__ == '__main__':
    d = read_m("軌道檔案.txt")
    tm = d["time"]
    x = d["x"]
    y = d["y"]
    z = d["z"]


    print("========迴歸係數=========")
    a = huiguixishu(tm, x)
    b = huiguixishu(tm, y)
    c = huiguixishu(tm, z)

    print(a)
    print(b)
    print(c)

    print("========預測=========")
    guji_time = [4200,4500,4800]
    beta0_list = [a[0],b[0],c[0]]
    beta1_list = [a[1],b[1],c[1]]

    for i in range(guji_time.__len__()):
        x = predic(guji_time[i],beta0_list[0],beta1_list[0])
        y = predic(guji_time[i],beta0_list[1],beta1_list[1])
        z = predic(guji_time[i],beta0_list[2],beta1_list[2])

        print(guji_time[i],format(x,'0.3f') ,format(y,'0.3f'),format(z,'0.3f'))

結果

========迴歸係數=========
[21146.959615384614, 1.2183738095238088]
[-7019.398461538461, 0.21143040293040288]
[15712.87576923077, -2.8401093406593407]
========預測=========
4200 26264.130 -6131.391 3784.417
4500 26629.642 -6067.962 2932.384
4800 26995.154 -6004.533 2080.351