数学建模python实现基础编程6

数学建模python实现基础编程6

1.monte carlo随机模拟1

a = 2.5
b = 5
u = 1
n = 10000
m1 = 0
for  i in range(1,400):#后面的数字根据感觉填
    d = np.random.normal(200,20,n)
    m2 = np.mean((b-a)*u*(u<=d) + ((b-a)*u-a*(u-d))*(u>d))
    if m2 > m1:
        m1 = m2
        u = u+1
    else:
        print("最佳值：",u-1)
        break

2.monte carlo随机模拟2

a = 2.5
b = 5
u = 1
n = 10000
m1 = 0
for  i in range(1,400):#后面的数字根据感觉填
    d = np.random.normal(200,20,n)
    m2 = np.mean((b-a)*u*(u<=d) + ((b-a)*u-a*(u-d))*(u>d))
    if m2 > m1:
        m1 = m2
        u = np.mean(d)
    else:
        print("甲方最佳值：",u)
        break
a = 2.5
b = 5
u = np.random.normal(200,20,n)
n = 10000
m1 = 0
d = np.random.normal(200,20,n)
aa = []
for  i in range(1,400):#后面的数字根据感觉填
    u = (d+u)/2
    m2 = np.mean((b-a)*u*(u<=d) + ((b-a)*u-a*(u-d))*(u>d))
    aa.append(np.mean(u))
    if m2 > m1:
        m1 = m2
        t = u
        u = d
        d = t
    else:
        print("乙方最佳值：",aa[-2])
        break

3.遗传算法解函数最大值

这个编程只能求函数值为正的

import numpy as np


def initpop(popsize, binarylength):
    pop = np.random.randint(0, 2, (popsize, binarylength))  # 生成popsize×binarylength的二维0、1序列
    return pop


def bintodec(ypop):
    pop = ypop.copy()
    [row, col] = pop.shape
    for i in range(col):
        pop[:, i] = 2 ** (col - 1 - i) * pop[:, i]
    pop = np.sum(pop, axis=1)
    num = []
    num = pop * 10 / 1023
    return num


# 计算种群适应度
def cal_objval(pop):
    x = bintodec(pop)
    objval = np.abs(x**3*np.cos(x))
    return objval


def selection(pop, fitval, popsize):
    idx = np.random.choice(np.arange(popsize), size=popsize, replace=True, p=fitval / fitval.sum())
    return pop[idx]


def crossover(pop, pc):
    [px, py] = pop.shape
    newpop = np.ones((px, py))
    for i in range(0, px, 2):
        if np.random.rand() < pc:
            cpoint = int(np.random.rand() * py * 10 // 10)
            newpop[i, 0:cpoint] = pop[i, 0:cpoint]
            newpop[i, cpoint:py] = pop[i + 1, cpoint:py]
            newpop[i + 1, 0:cpoint] = pop[i + 1, 0:cpoint]
            newpop[i + 1, cpoint:py] = pop[i, cpoint:py]
        #             newpop[i+1,:]=[pop[i+1,0:cpoint],pop[i,cpoint:py]]
        else:
            newpop[i, :] = pop[i, :]
            newpop[i + 1, :] = pop[i + 1, :]
    return newpop


def mutation(pop, pm):
    [px, py] = pop.shape
    newpop = np.ones((px, py))
    for i in range(px):
        if (np.random.rand() < pm):
            mpoint = int(np.random.rand() * py * 10 // 10)
            if mpoint <= 0:
                mpoint = 1
            newpop[i, :] = pop[i, :]
            if newpop[i, mpoint] == 0:
                newpop[i, mpoint] = 1
            else:
                newpop[i, mpoint] = 0
        else:
            newpop[i, :] = pop[i, :]
    return newpop


def best(pop, fitvalue):
    [px, py] = pop.shape
    bestindividual = pop[0, :]
    bestfit = fitvalue[0]
    for i in range(1, px):
        if fitvalue[i] > bestfit:
            bestindividual = pop[i, :]
            bestfit = fitvalue[i]
    return bestindividual, bestfit


if __name__ == "__main__":
    popsize = 100  # 种群规模
    binarylength = 10  # 二进制编码长度（DNA）
    pc = 0.6  # 交叉概率
    pm = 0.001  # 变异概率
    pop = initpop(popsize, binarylength)  # 初始化种群

    # 进行计算
    for i in range(100):
        # 计算当前种群适应度
        objval = cal_objval(pop)
        fitval = objval

        # 选择操作
        newpop = selection(pop, fitval, popsize)
        # 交叉操作
        newpop = crossover(newpop, pc)
        # 变异操作
        newpop = mutation(newpop, pm)
        # 更新种群
        pop = newpop

        # 寻找最优解并绘图
        [bestindividual, bestfit] = best(pop, fitval)

        x1 = bintodec(newpop)
        y1 = cal_objval(newpop)
        x = np.arange(-1.57, 20.18, 0.1)
        y = x**3*np.cos(x)
        if i % 10 == 0:
            plt.figure()
            plt.rcParams["font.sans-serif"] = ["SimHei"]  # 解决中文乱码问题
            plt.rcParams["axes.unicode_minus"] = False  # 使一些符号正常显示
            plt.plot(x, y)
            plt.plot(x1, y1, '*')
            plt.title('迭代次数为:%d' % i)
    [n] = bestindividual.shape
    x2 = 0
    for i in range(n):
        x2 += 2 ** (n - 1 - i) * bestindividual[i]
    print("The best X is ", x2 * 10 / 1023)
    print("The best Y is ", bestfit)

4.调用sklearn中的BP神经api使用预测

注意若是一元，则要把x,y内的数据分别括起来

import numpy as np
from sklearn.neural_network import MLPRegressor
x0 = np.array([[0],[34],[67],[101],[135],[202],[259],[336],[404],[471]])
y0 = np.array([[15.10],[21.36],[25.72],[32.29],[34.03],[39.45],[43.15],[43.46],[40.83],[30.75]])#由于是一元
md = MLPRegressor(solver='lbfgs',alpha=1e-5,hidden_layer_sizes=10)
md.fit(x0,y0)
print(md.score(x0,y0))
print(md.predict(np.array([[100]])))

5.AR模型

这个模型没太用明白

from statsmodels.tsa.ar_model import AR
d = np.array([10.010,11.260,9.000,9.090,9.440,9.090,8.730,8.680,9.040])
md = AR(d)
md_fit = md.fit()
print(md_fit.params)
print(md_fit.k_ar)
pred = md_fit.predict(start=7,end=len(d)+9,dynamic=False)
print(pred)

6.svm回归模型

import sklearn.svm as sv
x = np.array([50,60,60,70,70,80,80,90,90,90,95,100,100,100,105,105,110,110,110,115,115,115,120,120,120,125,130,130,135,135,140,140,145,150,150,155,155,160,160,160,165,170,180])
y = np.array([19,20,21,17,22,25,23,21,25,31,25,30,29,33,35,32,30,28,30,31,36,30,36,25,28,28,31,32,34,25,26,33,31,36,33,41,33,40,30,37,32,35,38])
x = x.reshape(-1,1)
md = sv.SVR(gamma='auto')
md.fit(x,y)
pred_y = md.predict([[120]])
print(pred_y)
score = md.score(x,y)
print(score)

效果还是不是很理想，误差还是有点大