Monte Carlo 使用 numpy 模拟三角和正态概率密度分布
Monte Carlo Simulation with triangular and normal probability density distibution using numpy
我正在尝试执行 Monte Carlo 模拟,以计算热泵系统电力成本的不确定性。我有几个输入参数(COP、电费),它们是三角概率分布。总电力成本由三个子组件(热泵和泵)的计算成本之和组成,并且服从(近似)正态概率分布。
我想知道我是否正确执行了 MC 模拟。由于我必须在 70 种不同的热泵系统上循环这个 MC 模拟,我也想知道是否有更快的方法。
由于我是编码方面的绝对新手,请为我的乱码道歉。
感谢您的帮助!
我的代码:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import triangular
N = 1_000_000
def energy_output(coef_performance, energy_input):
return energy_input * coef_performance / (coef_performance - 1)
COP_DISTRIBUTION_PARAM = dict(left=4, mode=4.5, right=5)
def seed_cop():
return triangular(**COP_DISTRIBUTION_PARAM )
INPUT_ENERGY_HEATING = 866
INPUT_ENERGY_COOLING = 912
def random_energy_output():
return energy_output(seed_cop(), energy_input=INPUT_ENERGY_HEATING)
energy_outputs = [random_energy_output() for _ in range(N)]
a = min(energy_outputs)
b = max(energy_outputs)
med = np.median(energy_outputs)
############################
def elec_costs_heatpump(elec_costs, coef_performance,energy_output):
return energy_output * 1000 / coef_performance * elec_costs
ELEC_DISTRIBUTION_PARAM = dict(left=0.14, mode=0.15, right=0.16)
def seed_elec():
return triangular(**ELEC_DISTRIBUTION_PARAM )
HP_OUTPUT_DISTRIBUTION_PARAM = dict(left=a, mode=med, right=b)
def seed_output():
return triangular(**HP_OUTPUT_DISTRIBUTION_PARAM )
def random_elec_costs_heatpump():
return elec_costs_heatpump(seed_elec(),seed_cop(),seed_output() )
elec_costs_heatpump = [random_elec_costs_heatpump() for _ in range(N)]
mean_hp = np.mean(elec_costs_heatpump)
std_hp = np.std(elec_costs_heatpump)
############################
def elec_costs_coldpump(elec_costs, coef_performance_pump,energy_input):
return energy_input * 1000 / coef_performance_pump * elec_costs
COP_PUMP_DISTRIBUTION_PARAM = dict(left=35, mode=40, right=45)
def seed_cop_pump():
return triangular(**COP_PUMP_DISTRIBUTION_PARAM )
def random_elec_costs_coldpump():
return elec_costs_coldpump(seed_elec(),seed_cop_pump(), energy_input=INPUT_ENERGY_COOLING)
elec_costs_coldpump = [random_elec_costs_coldpump() for _ in range(N)]
mean_cp = np.mean(elec_costs_coldpump)
sdt_cp = np.std(elec_costs_coldpump)
#########################
def elec_costs_warmpump(elec_costs, coef_performance_pump,energy_input):
return energy_input * 1000 / coef_performance_pump * elec_costs
def random_elec_costs_warmpump():
return elec_costs_warmpump(seed_elec(),seed_cop_pump(), energy_input=INPUT_ENERGY_HEATING)
elec_costs_warmpump = [random_elec_costs_warmpump() for _ in range(N)]
mean_wp = np.mean(elec_costs_warmpump)
sdt_wp = np.std(elec_costs_warmpump)
#########################
def total_costs(costs_heatpump, costs_coldpump, costs_warmpump):
return costs_heatpump + costs_coldpump + costs_warmpump
ELEC_COSTS_HEATPUMP_PARAM = dict(loc=mean_hp, scale=sdt_hp)
def seed_costs_hp():
return np.random.normal(**ELEC_COSTS_HEATPUMP_PARAM )
ELEC_COSTS_COLDPUMP_PARAM = dict(loc=mean_cp, scale=sdt_cp)
def seed_costs_cp():
return np.random.normal(**ELEC_COSTS_COLDPUMP_PARAM )
ELEC_COSTS_WARMPUMP_PARAM = dict(loc=mean_wp,scale=sdt_wp)
def seed_cost_wp():
return np.random.normal(**ELEC_COSTS_WARMPUMP_PARAM )
def random_total_costs():
return seed_costs_hp(), seed_costs_cp(), seed_cost_wp()
total_costs = [random_total_costs() for _ in range(N)]
print(total_costs)
#Plot = plt.hist(total_costs, bins=75, density=True)
太好了,你有一个代码原型!
对代码结构和可读性的一些印象:
最快的改进是将函数和脚本部分分开,这允许将您的算法拆分为简单的可测试块并将计算控制在一个地方,然后进行绘图
一些重复的代码可以转到自己的函数
坚持更广泛接受的命名约定 (PEP8) 是有回报的,这样人们就不会对风格感到惊讶,并且可以将更多注意力放在代码实质上。具体来说,您通常将函数命名为小写下划线 def do_something():
并且 UPPERCASE
是为常量保留的。
需要能够 运行 您的代码来检查加速,请参阅上面关于缺少什么的评论。
有关问题的更完整代码的更新,一些附加评论:
- 让函数成为函数,不要将函数参数作为全局变量传递
- 考虑将您的模型结构(您的方程式)与参数(固定输入和种子值)分开,它会在更长的时间内获得回报 运行
- 避免'magic numbers',将它们放入常量(例如
866
)或作为参数显式传递。
这里有一个需要考虑的重构:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import triangular
N = 1_000_000
#def EnergyHP():
# COP_Hp = triangular(4.0, 4.5, 5) # COP of heatpump
# Ein = 866 # Energy input heatpump
# return Ein * COP_Hp /(COP_Hp-1) # Calculation of Energy output of heatpump
#
#Eout = np.zeros(N)
#for i in range(N):
# Eout[i] = EnergyHP()
#minval = np.min(Eout[np.nonzero(Eout)])
#maxval = np.max(Eout[np.nonzero(Eout)])
#Medi= np.median(Eout, axis=0)
def energy_output(coef_performance, energy_input):
"""Return energy output of heatpump given efficiency and input.
Args:
coef_performance - <description, units of measurement>
energy_input - <description, units of measurement>
"""
return energy_input * coef_performance / (coef_performance - 1)
# you will use this variable again, so we put it into a function to recylce
COP_DISTRIBUTION_PARAM = dict(left=4, mode=4.5, right=5)
def seed_cop():
"""Get random value for coef of performance."""
return triangular(**COP_DISTRIBUTION_PARAM )
# rename it if it is a constant, pass as an argument is it is not
SOME_MEANINGFUL_CONSTANT = 866
def random_energy_output():
return energy_output(seed_cop(), energy_input=SOME_MEANINGFUL_CONSTANT)
# pure python list and metrics
energy_outputs = [random_energy_output() for _ in range(N)]
a = min(energy_outputs)
b = max(energy_outputs)
med = np.median(energy_outputs)
# above does this does not use numpy, but you can convert it to vector
energy_outputs_np = np.array(energy_outputs)
# or you can construct np array directly, this is a very readable way
# make a vector or random arguments
cop_seeded = triangular(**COP_DISTRIBUTION_PARAM, size=N)
# apply function
energy_outputs_np_2 = energy_output(cop_seeded, SOME_MEANINGFUL_CONSTANT)
下一个紧迫的意图是写一个 HeatPump
class,但我建议坚持
尽可能长时间地发挥作用——这通常会让我们想到 class
状态和方法更好。
我还发现参数的种子值可能不是独立的,在某些实验中您可以从联合分布中提取它们。
我正在尝试执行 Monte Carlo 模拟,以计算热泵系统电力成本的不确定性。我有几个输入参数(COP、电费),它们是三角概率分布。总电力成本由三个子组件(热泵和泵)的计算成本之和组成,并且服从(近似)正态概率分布。
我想知道我是否正确执行了 MC 模拟。由于我必须在 70 种不同的热泵系统上循环这个 MC 模拟,我也想知道是否有更快的方法。
由于我是编码方面的绝对新手,请为我的乱码道歉。
感谢您的帮助!
我的代码:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import triangular
N = 1_000_000
def energy_output(coef_performance, energy_input):
return energy_input * coef_performance / (coef_performance - 1)
COP_DISTRIBUTION_PARAM = dict(left=4, mode=4.5, right=5)
def seed_cop():
return triangular(**COP_DISTRIBUTION_PARAM )
INPUT_ENERGY_HEATING = 866
INPUT_ENERGY_COOLING = 912
def random_energy_output():
return energy_output(seed_cop(), energy_input=INPUT_ENERGY_HEATING)
energy_outputs = [random_energy_output() for _ in range(N)]
a = min(energy_outputs)
b = max(energy_outputs)
med = np.median(energy_outputs)
############################
def elec_costs_heatpump(elec_costs, coef_performance,energy_output):
return energy_output * 1000 / coef_performance * elec_costs
ELEC_DISTRIBUTION_PARAM = dict(left=0.14, mode=0.15, right=0.16)
def seed_elec():
return triangular(**ELEC_DISTRIBUTION_PARAM )
HP_OUTPUT_DISTRIBUTION_PARAM = dict(left=a, mode=med, right=b)
def seed_output():
return triangular(**HP_OUTPUT_DISTRIBUTION_PARAM )
def random_elec_costs_heatpump():
return elec_costs_heatpump(seed_elec(),seed_cop(),seed_output() )
elec_costs_heatpump = [random_elec_costs_heatpump() for _ in range(N)]
mean_hp = np.mean(elec_costs_heatpump)
std_hp = np.std(elec_costs_heatpump)
############################
def elec_costs_coldpump(elec_costs, coef_performance_pump,energy_input):
return energy_input * 1000 / coef_performance_pump * elec_costs
COP_PUMP_DISTRIBUTION_PARAM = dict(left=35, mode=40, right=45)
def seed_cop_pump():
return triangular(**COP_PUMP_DISTRIBUTION_PARAM )
def random_elec_costs_coldpump():
return elec_costs_coldpump(seed_elec(),seed_cop_pump(), energy_input=INPUT_ENERGY_COOLING)
elec_costs_coldpump = [random_elec_costs_coldpump() for _ in range(N)]
mean_cp = np.mean(elec_costs_coldpump)
sdt_cp = np.std(elec_costs_coldpump)
#########################
def elec_costs_warmpump(elec_costs, coef_performance_pump,energy_input):
return energy_input * 1000 / coef_performance_pump * elec_costs
def random_elec_costs_warmpump():
return elec_costs_warmpump(seed_elec(),seed_cop_pump(), energy_input=INPUT_ENERGY_HEATING)
elec_costs_warmpump = [random_elec_costs_warmpump() for _ in range(N)]
mean_wp = np.mean(elec_costs_warmpump)
sdt_wp = np.std(elec_costs_warmpump)
#########################
def total_costs(costs_heatpump, costs_coldpump, costs_warmpump):
return costs_heatpump + costs_coldpump + costs_warmpump
ELEC_COSTS_HEATPUMP_PARAM = dict(loc=mean_hp, scale=sdt_hp)
def seed_costs_hp():
return np.random.normal(**ELEC_COSTS_HEATPUMP_PARAM )
ELEC_COSTS_COLDPUMP_PARAM = dict(loc=mean_cp, scale=sdt_cp)
def seed_costs_cp():
return np.random.normal(**ELEC_COSTS_COLDPUMP_PARAM )
ELEC_COSTS_WARMPUMP_PARAM = dict(loc=mean_wp,scale=sdt_wp)
def seed_cost_wp():
return np.random.normal(**ELEC_COSTS_WARMPUMP_PARAM )
def random_total_costs():
return seed_costs_hp(), seed_costs_cp(), seed_cost_wp()
total_costs = [random_total_costs() for _ in range(N)]
print(total_costs)
#Plot = plt.hist(total_costs, bins=75, density=True)
太好了,你有一个代码原型!
对代码结构和可读性的一些印象:
最快的改进是将函数和脚本部分分开,这允许将您的算法拆分为简单的可测试块并将计算控制在一个地方,然后进行绘图
一些重复的代码可以转到自己的函数
坚持更广泛接受的命名约定 (PEP8) 是有回报的,这样人们就不会对风格感到惊讶,并且可以将更多注意力放在代码实质上。具体来说,您通常将函数命名为小写下划线
def do_something():
并且UPPERCASE
是为常量保留的。
需要能够 运行 您的代码来检查加速,请参阅上面关于缺少什么的评论。
有关问题的更完整代码的更新,一些附加评论:
- 让函数成为函数,不要将函数参数作为全局变量传递
- 考虑将您的模型结构(您的方程式)与参数(固定输入和种子值)分开,它会在更长的时间内获得回报 运行
- 避免'magic numbers',将它们放入常量(例如
866
)或作为参数显式传递。
这里有一个需要考虑的重构:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import triangular
N = 1_000_000
#def EnergyHP():
# COP_Hp = triangular(4.0, 4.5, 5) # COP of heatpump
# Ein = 866 # Energy input heatpump
# return Ein * COP_Hp /(COP_Hp-1) # Calculation of Energy output of heatpump
#
#Eout = np.zeros(N)
#for i in range(N):
# Eout[i] = EnergyHP()
#minval = np.min(Eout[np.nonzero(Eout)])
#maxval = np.max(Eout[np.nonzero(Eout)])
#Medi= np.median(Eout, axis=0)
def energy_output(coef_performance, energy_input):
"""Return energy output of heatpump given efficiency and input.
Args:
coef_performance - <description, units of measurement>
energy_input - <description, units of measurement>
"""
return energy_input * coef_performance / (coef_performance - 1)
# you will use this variable again, so we put it into a function to recylce
COP_DISTRIBUTION_PARAM = dict(left=4, mode=4.5, right=5)
def seed_cop():
"""Get random value for coef of performance."""
return triangular(**COP_DISTRIBUTION_PARAM )
# rename it if it is a constant, pass as an argument is it is not
SOME_MEANINGFUL_CONSTANT = 866
def random_energy_output():
return energy_output(seed_cop(), energy_input=SOME_MEANINGFUL_CONSTANT)
# pure python list and metrics
energy_outputs = [random_energy_output() for _ in range(N)]
a = min(energy_outputs)
b = max(energy_outputs)
med = np.median(energy_outputs)
# above does this does not use numpy, but you can convert it to vector
energy_outputs_np = np.array(energy_outputs)
# or you can construct np array directly, this is a very readable way
# make a vector or random arguments
cop_seeded = triangular(**COP_DISTRIBUTION_PARAM, size=N)
# apply function
energy_outputs_np_2 = energy_output(cop_seeded, SOME_MEANINGFUL_CONSTANT)
下一个紧迫的意图是写一个 HeatPump
class,但我建议坚持
尽可能长时间地发挥作用——这通常会让我们想到 class
状态和方法更好。
我还发现参数的种子值可能不是独立的,在某些实验中您可以从联合分布中提取它们。