Python中概率累计分布函数（CDF）分析

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

# 均值10，方差1,正态分布模拟数据
data = np.random.normal(10, 1, 100)

#计算正态概率密度函数在x处的值
def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=np.mean(data), scale=np.std(data))
    return y

#计算正态分布累积概率值
def norm_dist_cdf(theta):
    y = norm.cdf(theta,loc=10,scale=1)
    return y
    
#利用ppf找到适合的横坐标，百分点函数
#ppf分位点函数（CDF的逆）即累计分布函数的逆函数（分位点函数，给出分位点返回对应的x值）。
#scipy.stats.norm.ppf(0.95， loc=0，scale=1)返回累积分布函数中概率等于0.95对应的x值（CDF函数中已知y求对应的x）。
x = np.linspace(norm.ppf(0.01,loc=np.mean(data), scale=np.std(data)),
                norm.ppf(0.99,loc=np.mean(data), scale=np.std(data)), len(data))  #linspace() 函数返回指定间隔内均匀间隔数字的 ndarray。

y1=norm_dist_prob(x)
y2=norm_dist_cdf(x)
plt.plot(x, y1,'g',label='pdf')
plt.plot(x, y2,'r',label='cdf')
plt.show()

#CDF数据处理
def Cumsum_cdf(DATA):
    denominator = len(DATA['VALS'])
    Data1 = pd.Series(DATA['VALS'])

    # #利用value_counts方法进行分组频数计算
    Fre=Data1.value_counts()
    # #对获得的表格整体按照索引自小到大进行排序
    Fre_sort=Fre.sort_index(axis=0,ascending=True)
    # # 每个数据出现频数除以数据总数才能获得该数据的概率
    # #重置表格索引
    Fre_df=Fre_sort.reset_index()
    # #将频数转换成概率
    Fre_df['VALS']=Fre_df['VALS']/denominator
    # #将列表列索引重命名
    Fre_df.columns=['Rds','Fre']

    # # 将数据列表从小到大排列，然后将每个数据出现的概率进行叠加
    # #利用cumsum函数进行概率的累加并按照顺序添加到表格中
    Fre_df['cumsum']=np.cumsum(Fre_df['Fre'])

    return Fre_df
    
def Cumulative_Distribution_Function(DATA,VAL1 =300):
    import matplotlib.pyplot as plt
    from scipy.special import jn, jn_zeros
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

    df_A1,df_A2,df_A3 = DATA
    #cdf数据处理函数,添加三条曲线
    Fre_A1 = Cumsum_cdf(df_A1[df_A1.VALS > 0])
    Fre_A2 = Cumsum_cdf(df_A2[df_A2.VALS > 0])
    Fre_A3 = Cumsum_cdf(df_A3[df_A3.VALS > 0])

    # 创建画布，只有一张图，也可以多张
    # plot = plt.figure()
    fig, ax1= plt.subplots(figsize=(8, 4), dpi=200)
    ax1.spines["right"].set_visible(False)
    ax1.spines["top"].set_visible(False)
    ax1.grid(ls="--", lw=0.25, color="#4E616C")

    # 按照Rds列为横坐标，累计概率分布为纵坐标作图
    #添加三条数据曲线
    ax1.plot(Fre_A1['Rds'], Fre_A1['cumsum'], label=f'工况1', mfc="white", ms=5)
    ax1.plot(Fre_A2['Rds'], Fre_A2['cumsum'], label=f'工况2', mfc="white", ms=5)
    ax1.plot(Fre_A3['Rds'], Fre_A3['cumsum'], label=f'工况3', mfc="white", ms=5)


    # # 图的标题
    ax1.legend(loc='upper left',
              frameon=False,
              edgecolor="None",fontsize=8)

    #峰值负载
    upper_peak = max(Fre_A1['Rds'].max(),Fre_A2['Rds'].max(),Fre_A3['Rds'].max())

    f1_A25 = Fre_A3[Fre_A3['Rds'] / VAL1 >= 0.25]["cumsum"].to_list()[0]
    f1_A50 = Fre_A3[Fre_A3['Rds'] / VAL1 >= 0.5]["cumsum"].to_list()[0]

    # # 显示坐标点横线、竖线
    peak = round(upper_peak/VAL1,2)*100
    
    # # 峰值点负载率
    ax1.annotate(
        f"$({peak}$%$,100$%$)$", # 标注文字
        (upper_peak,1),
        # size="small",#medium
        xytext=(-65, -40),
        textcoords="offset points",
        fontsize=10,
        arrowprops=dict(
            arrowstyle="->",
            connectionstyle="arc3,rad=-0.3"), )

    #25%、50%分界线
    plt.vlines(VAL1 * 0.25, 0,f1_A25, colors="c", linestyles="dashed")

    plt.vlines(VAL1 * 0.5, 0, f1_A50, colors="c", linestyles="dashed")

    #峰值线
    plt.vlines(upper_peak, 0,1,colors="r", linestyles="dashed")


    #添加阴影区域，工况1
    X = Fre_A1['Rds'].tolist()
    C = Fre_A1['cumsum'].tolist()
    x = np.array([i for i in X])
    plt.fill_between(x,C,where=(VAL1*0.5*0.8 < x) & (x< VAL1*0.5),facecolor='green')

    #标注工况3,50%负载概率值
    ax1.annotate(
        f"$(50$%$,{round(f1_A50,2)*100}$%$)$", # 标注文字
        (VAL1*0.5,f1_A50),
        # size="small",#medium
        xytext=(-65, -40),
        textcoords="offset points",
        fontsize=10,
        arrowprops=dict(
            arrowstyle="->",
            connectionstyle="arc3,rad=-0.3"), )

    #标注分界线名称
    ax1.text(
        VAL1 * 0.27, 0.10,
        f'25%负载率界限',
        va="top", ha="left",
        size=8)

    ax1.text(
        VAL1 * 0.52, 0.10,
        f'50%负载率界限',
        va="top", ha="left",
        size=8)

    #主题名称
    team_ = ax1.text(
        x=0, y=ax1.get_ylim()[1] + ax1.get_ylim()[1] / 20,
        s=f"Cumulative probability",
        color="#4E616C",
        va='center',
        ha='left',
        size=8
    )
    # 横轴名
    ax1.set_xlabel("用电量（kWh）",fontsize =8)

    plt.savefig(r".\test.png",bbox_inches="tight")

    return plt

# 外部导入数据
DF = pd.read_excel(r".\ST\分析数据.xlsx")
df_A1 = DF[["工况1"]].copy()
df_A1 = df_A1.rename(columns={'工况1': 'VALS'})
df_A2 = DF[["工况2"]].copy()
df_A2 = df_A2.rename(columns={'工况2': 'VALS'})
df_A3 = DF[["工况3"]].copy()
df_A3 = df_A3.rename(columns={'工况3': 'VALS'})
DATA = df_A1,df_A2,df_A3