数据统计与分析基础实验（三）Hebuter慎抄

实验三 常规数学统计运算

实验目的：掌握各种常见的数学统计、数学分析运算的编程实现过程。

实验内容：

1、随机生成一个10x15的高斯矩阵，均值为自己学号后两位，方差为1。对该矩阵分别进行LU、QR、奇异值，并展示分解结果。

2、利用软件求解优化问题：

min Z = - X1 - 0.8X2 - 1.2X3, 

约束：X1 –X2 + X3 <= 30,

3X1 + 2X2 + 4X3 <=42,

3X1 + 2X2 <= 30,

X1, X2 , X3 >=0。

    并找到其中有效约束。

3、对于函数f(x) = (x2-5x+3)exp(-4x)cosx, x = 0: 0.1:1,进行三种方式插值，并将插值曲线与原曲线绘制在同一个Figure窗口。

4、统计自己过去12个月实际生活花费的数值，并拟合成一条一次、二次、三次曲线，三条曲线分别用不同颜色和线型展示在同一张图里。

5、求解三重积分∫（-0.5，1）∫（0，0.5）∫（-0.5pi,pi）(y sinx exp(x)+z cosx x2)dxdydz.

6、求微分方程 xy’+ y – exp(x) = 0,在初始条件y(1) = 2 exp下的特解，并画出解函数的图形。

7、 某人进行射击，及每次命中的概率为0.1，独立射击50次，求击中10次以上且40次以下的概率。

8、假设自己过去12个月实际生活花费的数值服从正态分布，请求出其均值和方差的极大似然估计。

实验代码及结果：

1、

import numpy as np
import sympy
import xlwt

matrix = np.random.normal(loc=98, scale=1, size=(10, 15))


def save(array, name):
    filename = xlwt.Workbook()  # 创建工作簿
    sheet = filename.add_sheet(u'sheet1', cell_overwrite_ok=True)  # 创建sheet

    [h, l] = array.shape  # h为行数，l为列数
    for i in range(h):
        for j in range(l):
            sheet.col(j).width = 3300
            sheet.write(i, j, float(array[i, j]))
    filename.save(name)


def LU(m):
    l, u, _ = m.LUdecomposition()
    print(np.array(l))
    print(np.array(u))
    save(m, 'matrix.xls')
    save(l, 'l_matrix.xls')
    save(u, 'u_matrix.xls')


def QR(m):
    Q, R = m.QRdecomposition()
    print(np.array(Q))
    print(np.array(R))
    save(Q, 'Q_matrix.xls')
    save(R, 'R_matrix.xls')


def SVD(m):
    S = m.singular_value_decomposition()
    for i in range(len(S)):
        print(np.array(S[i]))
        save(S[i], 'SVD_matrix_' + str(i) + '.xls')


def main():
    m = sympy.Matrix(matrix)
    print(np.array(m))
    LU(m)
    QR(m)
    SVD(m)


if __name__ == '__main__':
    main()

 2、

from scipy.optimize import minimize
import numpy as np


def fun(x):
    return - (x[0] + 0.8 * x[1] + 1.2 * x[2])


if __name__ == '__main__':
    # 约束
    cons = ({'type': 'ineq', 'fun': lambda x: 30 - x[0] + x[1] - x[2]},
            {'type': 'ineq', 'fun': lambda x: 42 - 3 * x[0] - 2 * x[1] - 4 * x[2]},
            {'type': 'ineq', 'fun': lambda x: 30 - 3 * x[0] - 2 * x[1]},
            {'type': 'ineq', 'fun': lambda x: x[0]},
            {'type': 'ineq', 'fun': lambda x: x[1]},
            {'type': 'ineq', 'fun': lambda x: x[2]}
            )
    x0 = np.array((0, 0, 0))  # 设置初始值
    res = minimize(fun, x0, options={'disp': True}, constraints=cons)
    print(res)
    print('最小值：', res.fun)
    print('最优解：', res.x)
    print('迭代终止是否成功：', res.success)
    print('迭代终止原因：', res.message)

3、

import numpy as np
import math
from matplotlib import pyplot as plt
from scipy import interpolate
import pylab as pl

x = np.arange(0, 1.1, 0.1)


def fun(x):
    y = []
    for i in range(len(x)):
        y.append((x[i] * x[i] - 5 * x[i] + 3) * math.exp(-4 * x[i]) * np.cos(x[i]))
    return np.array(y)


if __name__ == '__main__':
    print(x)

    xnew = np.arange(0, 1, 0.01)
    print(xnew)
    for kind in ["nearest", "slinear", "quadratic"]:  # 插值方式
        # "nearest"为阶梯插值
        # slinear 线性插值
        # "quadratic"为2阶样条曲线插值
        f = interpolate.interp1d(x, fun(x), kind=kind)
        # ‘slinear', ‘quadratic' and ‘cubic' refer to a spline interpolation of first, second or third order)
        ynew = f(xnew)
        pl.plot(xnew, ynew, label=str(kind))
    pl.plot(xnew, fun(xnew), linestyle="dashdot", color="black",label="origin")
    pl.legend(loc="lower right")
    pl.show()

 

 4、

import xlrd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from numpy import polyfit, poly1d


def read_data(filepath):
    table = xlrd.open_workbook(filepath).sheets()[0]
    data = table.col_values(0)
    print(data)
    X = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
    return X, data


def fit(X, Y):
    w = polyfit(X, Y, 1)
    print("Parameters of polynomial:")
    print(w)
    y = w[0] * X + w[1]
    plt.plot(X, y,label="first order")

    w = polyfit(X, Y, 2)
    print("Parameters of polynomial:")
    print(w)
    y = w[0] * X ** 2 + w[1] * X + w[2]
    plt.plot(X, y,label="second order")

    w = polyfit(X, Y, 3)
    print("Parameters of polynomial:")
    print(w)
    y = w[0] * X ** 3 + w[1] * X ** 2 + w[2] * X + w[3]
    plt.plot(X, y,label="third order")


if __name__ == '__main__':
    X, Y = read_data("E:\数据分析与统计基础\实验\Project\exp3\expend.xls")
    X = np.array(X)
    plt.scatter(X, Y)
    fit(X, Y)
    plt.legend(loc=2)
    plt.show()

 5、

import math
from scipy import integrate
import numpy as np

if __name__ == '__main__':
    f = lambda x, y, z: y * np.sin(math.exp(x)) + z * np.cos(x) * x * x
    v, err = integrate.tplquad(f, -0.5, 1, 0, 0.5, -0.5 * np.pi, np.pi)
    print("result:"+str(v))

6、

import math
from sympy import *
import matplotlib.pyplot as plt
import numpy as np

if __name__ == '__main__':
    f = symbols('f', cls=Function)
    x = symbols('x')
    eq = Eq(x * f(x).diff(x, 1) + f(x) - exp(x), 0)
    print(dsolve(eq, f(x)))
    C1 = symbols('C1')
    eqr = (C1 + exp(x))/x
    # print(type(eqr))
    '''eqr = str(dsolve(eq,f(x))).split(', ')[1].strip(')')
    eqr.strip()
    print('s'+eqr)
    print(eqr)'''
    eqr1 = eqr.subs(x, 1)
    print(solveset(eqr - 2*math.e, C1))
    eqr2 = eqr.subs(C1, 5.43656365691809 - E)
    print(eqr2)
# 画图

    x_1 = np.arange(0.001, 5, 0.1)
    x_2 = np.arange(-5,-0.001,0.1)
    y_1 = [(E + exp(x)/x) for x in x_1]
    y_2 = [(E + exp(x) / x) for x in x_2]
    plt.plot(x_1, y_1)
    plt.plot(x_2, y_2)
    plt.axis([-5,5,-10,10])
    plt.grid()
    plt.show()

 

 

7、

import math
from decimal import *

import numpy as np


# getcontext().prec = 100
def C(n, m):
    return math.factorial(n) / math.factorial(m) / math.factorial(n - m)


def main():
    P = []
    for i in range(11, 40):
        P.append(C(50, i) * 0.1 ** i * 0.9 ** (50 - i))
        # print(C(50, i) * 0.1 ** i * 0.9 ** (50 - i))
        # print(np.sum(np.array(P)))
    print("概率为:%.8f"%np.sum(np.array(P)))


if __name__ == '__main__':
    main()

 

8、

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
import xlrd


def read_data(filepath):
    table = xlrd.open_workbook(filepath).sheets()[0]
    data = table.col_values(0)
    print("原数据为:")
    print(data)
    return data


if __name__ == '__main__':
    Y = read_data("E:\数据分析与统计基础\实验\Project\exp3\expend.xls")
    print("均值方差分别为：")
    print(norm.fit(Y))