该篇介绍建模常用算法，主要使用sicpy

from scipy import optimize

使用sicpy.optimize模块的linprog函数可以求解线性规划问题，该函数中规划的标准型如下:

rse = scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='simplex', callback=None, options=None)

参数:

• c: 目标向量(目标函数系数)
• A_ub, b_ub: 不等号约束(默认为<=，可以取符号得到>=号)
• A_eq, b_eq: 等号约束
• bounds: 决策向量是上下界的n个元素组成的数组
• method: 求解器的类型，如单纯形法(更多请查看官方文档)
• callback: 调用回调函数，基本不使用也不懂

返回值:

• res.x: 最优解(不等式未知数的解)
• res.fun: 目标函数最优值(要求的最值)

也可以使用cvxopt.slovers模块求解线性规划问题

import numpy as np # 必要
from cvxopt import materix, slovers

该线性规划模型标准型:

使用方法类似:

c = materix([]) # 目标向量
A = materix([])
b = materix([]) # 不等号约束

# Aqe, beq等号约束，同理

sol = solvers.lp(c, A, b, Aeq, beq)

• 参数要用materix类型
• sol['x'] 为最优解
• sol['primal objective'] 为最优值

cvxpyt是一个用于凸优化的模块，凸优化为优化的一种类型，包括线性规划(LP)，二次规划(Quadratic Program)， 二次约束的二次规划(Quadratically Contrained Quaratic Program)，半正定规划(SDP, Semidefinite Program)等，关于凸优化具体可以看这篇博客，这是cvxpy官网
使用cvxpy求解凸优化问题

import cvxpy as cp

在CVXPY中变量有标量(只有数值大小)，向量，矩阵，也有常量(Parameter)
用例:

# Create two scalar optimization variables.
# 定义两个标量
x = cp.Variable()
y = cp.Variable()

# Create two constraints.
# 创建一个约束式
constraints = [x + y == 1,
               x - y >= 1]

# Form objective.
# 生成目标函数
obj = cp.Minimize((x - y)**2)

# Form and solve problem.
# 生成规划问题
prob = cp.Problem(obj, constraints)

# Returns the optimal value.
# 返回最优解
prob.solve() 

print("status:", prob.status)
print("optimal value:", prob.value)
print("optimal var:", x.value, y.value)

• prob.status 为求解状态
  optimal: 问题被成功解决
  infeasible：问题无解
  unbounded：无边界
  optimal_inaccurate：解不精确

  如果CVXPY求解器求解完全失败，CVXPY将会抛出一个SolverError异常，如果这发生，你应该使用其他cvxpy求解器

• prob.value 为目标函数最优值

• x.value, y.value 为最优解

说人话就是很多情况下x的系数是会变的，要分析这些变化对最优解的影响

# 就是那个牛奶生产问题
c = [-72, -64]
A = [[1, 1], [12, 8]]
b = [[50], [480]]
bound = ((0, 100/3.0), (0, None))
res = optimize.linprog(c, A, b, None, None, bound, method='simplex', options={"disp": True})
print(res)

输出

Optimization terminated successfully.
Current function value: -3360.000000
Iterations: 4
     con: array([], dtype=float64)
     fun: -3360.0
 message: 'Optimization terminated successfully.'
     nit: 4
   slack: array([0., 0.])
  status: 0
 success: True
       x: array([20., 30.])

其中，slack为两个0，说明约束条件为紧约束，即最优解的约束达到了边界，所以当约束边界变大时，最优解会更小(相反更大)

更多分布函数见scipy.stats

标准正态分布密度函数 官网文档

from scipy.stats import norm

$$f(x) = \frac{exp(-x^2/2)}{\sqrt{2\pi}}$$

• 分布函数 (Cumulative Distribution Function)

  norm.cdf(x, loc=0, scale=1)

1. 设$X\sim N(3,5^2)$，求$P(2<X<6)$?

   norm.cdf(6, 3, 5) - norm.cdf(2, 3, 5)

2. $P(-3c&lt;X&lt;2c)=0.6$，求$c$?

   from scipy.optimize import fsolve
   f = lambda c: norm.cdf(2*c, 3, 5) - norm.cdf(-3*c, 3, 5) - 0.6
   fsolve(f, 0)

• 概率密度函数（Probability Density Function）

  norm.pdf(x, loc=0, scale=1)

• 百分点函数，概率密度函数的积分值 (Percent Point Function)

  norm.ppf(q, loc=0, scale=1)

求标准正态分布下，$\alpha = 0.1$时，上$\alpha$分位数?

norm.ppf(1-0.001, 0, 1) # 上分位数

画出结果

x = np.linspace(-4, 4, 100) # x范围
y = norm.pdf(x, 0, 1)       # 计算标准正态分布的y值
plt.plot(x, y)

x2 = np.linspace(norm.ppf(1 - 0.1, 0, 1), 4, 100)
y2 = norm.pdf(x2)
plt.fill_between(x2, 0, y2, facecolor='orange') # 填充
plt.plot([-4, 4], [0, 0], color='black')

plt.text(1.9, 0.07, '$\leftarrow \\alpha $=0.1', size=14) # 注释
plt.show()

• 计算(二项)分布的均值、方差、偏度和峰度

  from scipy.stats import binom
  binom.stats(n, p, moments = 'mvsk')

  Eg: 计算b(20,0.8)的均值、方差、偏度和峰度

  mean, variance, skewness, kurtosis = binom.stats(20, o.8, moments = 'mvsk')