基于Sigma调整正态分布方差异常：当前输出与预期不符

我需要为多轮实验生成正态分布随机数组并写入CSV文件，同时打印每轮的均值和标准差。设置的均值mu=50.0、标准差sigma=10.0，但当前输出的标准差（代码中标注为var）约为1.0，与预期的10.0不符。

我的代码 #

from scipy.stats import truncnorm
import numpy as np
import os
import csv 
import pandas as pd
import random
import matplotlib.pyplot as plt

Runs=5

for i in range(0,Runs): 
    
    mu, sigma = 50.0, 10.0 # mean and standard deviation
    Nodes=220
    r = (1e-6)*np.random.normal(mu, sigma, Nodes)

    sort_r = np.sort(r)
    r1=sort_r[::-1]
    r1=r1.reshape(1,Nodes)
    r2 = r.copy()
    np.random.shuffle(r2.ravel()[1:])
    r2=r2.reshape(1,Nodes)
    maximum = r2.max()
    indice1 = np.where(r2 == maximum)

    r2[indice1] = r2[0][0]
    r2[0][0] = maximum

    r2[0][Nodes-1] = maximum 

    a = (1e-6) * (50 + (r - r.mean()) / r.std())
    mean=a.mean()*1e6
    var=a.std()*1e6
    print("Mean =", mean)
    print("var =", var)
    
    with open(rf"Radius info_{Nodes}ND.txt", 'w+') as f: 

        f.write(f"mean = {str(mean)}\n")
        f.write(f"var = {str(var)}\n")
        #f.write(f"std = {str(std*1e6)}\n")

    r=pd.Series(r*1e6)
    #print("r =",r)
    r.hist()
    plt.xlabel("Pore throat size [\u03BCm]",size=15)
    plt.ylabel("Count",size=15)
    plt.title(rf"var={round(var,1)}",size=25)
    plt.savefig(rf"C:\Users\USER\OneDrive - Technion\Research_Technion\Python_PNM\Surfactant A-D\220 nodes_2\Histogram.png")
    
    
def function(run):
    parent_folder = str(run + 1)
    os.mkdir(parent_folder)

         
    with open(os.path.join(parent_folder, rf"Inv_Radius_{Nodes}_{round(mean,1)}_{round(var,1)}ND.csv"), 'w+') as f:
        np.random.shuffle(r2[0][1:-1])
          #print(f'{run}:{r2}')
        inv_r=1/r2
          #print("inv_r =",[inv_r])
        writer = csv.writer(f)
        writer.writerows(inv_r)
     
for x in range(Runs): 
    function(x)

当前输出 #

Mean = 49.99999999999999
var = 1.0
Mean = 49.99999999999999
var = 1.0
Mean = 49.99999999999999
var = 0.9999999999999996
Mean = 49.99999999999999
var = 1.0000000000000004
Mean = 49.99999999999999
var = 1.0

预期输出 #

Mean = 49.99999999999999
var = 10.0
Mean = 49.99999999999999
var = 10.0
Mean = 49.99999999999999
var = 0.9999999999999996
Mean = 49.99999999999999
var = 10.0000000000000004
Mean = 49.99999999999999
var = 10.0

问题原因 #

核心问题出在这段代码：

a = (1e-6) * (50 + (r - r.mean()) / r.std())

这里对原始数组r做了标准化处理——(r - r.mean())/r.std()会把数据转换成均值为0、标准差为1的分布，后续仅加上50再缩放，最终a的标准差自然是1左右（乘以1e6后显示为1.0），完全偏离了预期的10.0。

修复方案 #

方案1：直接计算原始数组的均值和标准差

你生成的r已经是符合mu=50、sigma=10的正态分布（转换为米单位），直接转换回微米计算即可：

# 替换原有的a、mean、var计算代码
mean = r.mean() * 1e6
var = r.std() * 1e6  # 注意：代码中你把标准差标注为var，若需方差应使用r.var() * 1e12

方案2：若需严格校准数组的均值和标准差

如果想修正原始样本的微小波动，确保数组严格符合目标均值和标准差，应调整标准化后的缩放逻辑：

a = (1e-6) * (mu + sigma * (r - r.mean()) / r.std())
mean = a.mean() * 1e6
var = a.std() * 1e6

这样标准化后乘以目标标准差sigma，再加上目标均值mu，就能得到符合预期的分布。

额外优化点 #

1. 文件覆盖问题：循环中每次会覆盖Radius info_{Nodes}ND.txt和直方图文件，建议文件名添加轮次标识，比如Radius info_{Nodes}ND_run{i+1}.txt、Histogram_run{i+1}.png
2. 函数变量错误：原function函数使用的mean和var是循环最后一轮的变量，需将这两个变量作为参数传入函数，避免闭包导致的错误
3. 绘图重叠问题：每轮绘图后需调用plt.close()关闭画布，防止多轮图像重叠

修复后的完整代码 #

from scipy.stats import truncnorm
import numpy as np
import os
import csv 
import pandas as pd
import random
import matplotlib.pyplot as plt

Runs=5
Nodes=220
mu, sigma = 50.0, 10.0 # mean and standard deviation

def function(run, r2, mean, var):
    parent_folder = str(run + 1)
    os.mkdir(parent_folder)
    with open(os.path.join(parent_folder, rf"Inv_Radius_{Nodes}_{round(mean,1)}_{round(var,1)}ND.csv"), 'w+') as f:
        np.random.shuffle(r2[0][1:-1])
        inv_r=1/r2
        writer = csv.writer(f)
        writer.writerows(inv_r)

for i in range(Runs): 
    # 生成符合要求的正态分布数组（米单位）
    r = (1e-6)*np.random.normal(mu, sigma, Nodes)

    sort_r = np.sort(r)
    r1=sort_r[::-1]
    r1=r1.reshape(1,Nodes)
    r2 = r.copy()
    np.random.shuffle(r2.ravel()[1:])
    r2=r2.reshape(1,Nodes)
    maximum = r2.max()
    indice1 = np.where(r2 == maximum)

    r2[indice1] = r2[0][0]
    r2[0][0] = maximum
    r2[0][Nodes-1] = maximum 

    # 直接计算原始数组的均值和标准差（转换回微米）
    mean = r.mean() * 1e6
    var = r.std() * 1e6
    print("Mean =", mean)
    print("var =", var)
    
    # 保存每轮信息，文件名添加轮次
    with open(rf"Radius info_{Nodes}ND_run{i+1}.txt", 'w+') as f: 
        f.write(f"mean = {str(mean)}\n")
        f.write(f"var = {str(var)}\n")

    # 绘制直方图，文件名添加轮次
    r_um = pd.Series(r*1e6)
    r_um.hist()
    plt.xlabel("Pore throat size [μm]",size=15)
    plt.ylabel("Count",size=15)
    plt.title(rf"var={round(var,1)}",size=25)
    plt.savefig(rf"C:\Users\USER\OneDrive - Technion\Research_Technion\Python_PNM\Surfactant A-D\220 nodes_2\Histogram_run{i+1}.png")
    plt.close()  # 关闭画布避免重叠
    
    # 传入当前轮次的变量到函数
    function(i, r2, mean, var)

内容的提问来源于stack exchange，提问作者KeplerNick123