Tech

物理实验工具开发小记

by illlights, 2023-03-24


初衷

物理实验有大量繁琐的数据处理,这些东西多而且重复,显然就是计算机该做的事,想要开发一个计算机程序理所当然。

其实学校已经有学长开发出针对不同实验的具体计算程序,但是缺乏计算过程、计算明显有问题,而且输入限制太大。为了能自己完整地掌控计算细节并且实现各种实验通用的部分,还是选择了自己重新做一个。

有AI的帮助,开发这种小东西基本就不在话下,鉴于AI比较“擅长”高层级一些的语言,所以我基本上是想都没想上 Python ,先告诉AI什么是 3σ 法去除坏值,然后一句“请用 Python 实现”,此时的我还没有发现自己开了一个并不简单的坑。

AI的成果

ChatGPT-4 整体来说做的也还不错,至少对于 3σ 去除坏值这种广泛使用的算法AI完成的还不错,第一遍生成的代码如下。

import numpy as np

def read_data(file_path):
    with open(file_path, 'r') as f:
        data_str = f.read()
    data_list = [float(x) for x in data_str.split()]
    return np.array(data_list)

def remove_outliers(data):
    mean = np.mean(data)
    std_dev = np.std(data)
    lower_bound = mean - 3 * std_dev
    upper_bound = mean + 3 * std_dev

    good_values = [x for x in data if lower_bound <= x <= upper_bound]
    outliers = [x for x in data if x not in good_values]

    return np.array(good_values), np.array(outliers)

def save_results(good_values, outliers, mean, std_dev):
    individual_std_devs = [abs(x - mean) for x in good_values]

    with open('result.txt', 'w') as f:
        f.write('剩余的值: ' + ' '.join(map(str, good_values)) + '\n')
        f.write('剔除的值: ' + ' '.join(map(str, outliers)) + '\n')
        f.write('最终的平均数和σ: ' + str(mean) + ' ' + str(std_dev) + '\n')
        f.write('每个good_values的标准差: ' + ' '.join(map(str, individual_std_devs)) + '\n')

if __name__ == '__main__':
    file_path = 'input.txt' # 修改为你的input.txt文件路径
    data = read_data(file_path)
    good_values, outliers = remove_outliers(data)
    mean = np.mean(good_values)
    std_dev = np.std(good_values)
    save_results(good_values, outliers, mean, std_dev)

但是它的缺点是,只会按照大多数人的写法写,基本上不能完美符合特殊化需求。比如这个程序使用的标准差直接调用 Numpy 库函数,std_dev = np.std(data),不仅少了置信系数,而且难以找出里面的详细数学处理,显然与我的初衷有些许偏差。这个问题的解决后面会详细说。

如果说这上面只是小问题,那计算不确定度的程序写的就可以说是完全错误了。

def calculate_uncertainty(data, instrument_error):
    mean = np.mean(data)
    std_dev = np.std(data)
    a_class_uncertainty = std_dev / np.sqrt(len(data))
    b_class_uncertainty = instrument_error
    combined_uncertainty = np.sqrt(a_class_uncertainty**2 + b_class_uncertainty**2)
    return mean, a_class_uncertainty, b_class_uncertainty, combined_uncertainty

首先还是公式存在问题,A 类不确定度和 B类不确定度都少了系数,在强调置信系数之后变成了这样

def calculate_uncertainty(data, instrument_error, confidence_coefficient):
    mean = np.mean(data)
    std_dev = np.std(data)
    a_class_uncertainty = std_dev / np.sqrt(len(data))
    b_class_uncertainty = instrument_error
    combined_uncertainty = np.sqrt(a_class_uncertainty**2 + b_class_uncertainty**2 / confidence_coefficient)
    return mean, a_class_uncertainty, b_class_uncertainty, combined_uncertainty

显然 AI 完全没有理解这些概念的实际意义,也没有理解这些公式到底是什么(这个问题在我问np.std()的具体计算过程时也有所体现),只是组合上去然后写出来。

可以这样说,至少在科学领域,AI 连最基本的工作都无法胜任。

完善

实际上只要花时间一点一点教 AI 改代码,它还是能把特殊化需求写好。但是,有这时间教 AI 写,为什么不自己改一下核心部分?我觉得这个问题是很多 AI 使用者所忽略的:既然 AI 只是一个工具,那么如果使用它花费的时间反而变多,应该果断舍弃这一低效的方法。

在查阅了大量的文档之后,发现 Numpy 本身就没有给出内置函数的数学表达式,只能查它的源代码——但是太慢了,作为一个本身误差就大的处理工具,基本正确就好。再查了一些资料之后,程序改为了

def remove_outliers(data):
    all_outliers = []
    
    while True:
        mean = np.mean(data)
        std_dev = np.std(data, ddof=1) * t683(data)
        lower_bound = mean - 3 * std_dev
        upper_bound = mean + 3 * std_dev

        good_values = [x for x in data if lower_bound <= x <= upper_bound]
        outliers = [x for x in data if x not in good_values]

        if len(outliers) == 0:
            break

        all_outliers.extend(outliers)
        data = np.array(good_values)

    return np.array(good_values), np.array(all_outliers)

这里的 while True:是为了让程序不断去除坏值直到无坏值,这种简单的逻辑 AI 也不能一次解决。还有t683()按书上给的值做好了对应。

不确定度的计算部分我基本上完全重写,本身也就是简单打个公式,只不过参数多一些参数的关系复杂很多。

在这之后,为了抄过程方便,我重写了输出格式,把两部分合成一个文件一次性输出。虽然这样写好像没有 AI 写的规范,但是这样一个简单的程序又没什么理解难度,怎么方便怎么来吧。

def save_results(good_values, outliers, mean, std_dev, instrument_error, confidence_coefficient):
    individual_std_devs = [abs(x - mean) for x in good_values]

    with open('result.txt', 'w') as f:
        f.write('剩余的值: ' + ' '.join(map(str, good_values)) + '\n')
        f.write('剔除的值: ' + ' '.join(map(str, outliers)) + '\n')
        f.write('最终的平均数: ' + str(mean) + '\n')
        f.write('S = ' + str(np.std(good_values, ddof=1)) + '\n')
        f.write('σ = '+ str(std_dev) + ', 3σ = '+ str(3*std_dev) + '\n')
        f.write('每个最终值的标准差: ' + ' '.join(map(str, individual_std_devs)) + '\n')
        a_class_uncertainty = std_dev / np.sqrt(np.size(good_values))
        b_class_uncertainty = instrument_error / confidence_coefficient
        combined_uncertainty = np.sqrt(a_class_uncertainty**2 + b_class_uncertainty**2)
        f.write('A类不确定度(标准偏差): ' + str(a_class_uncertainty) + '\n')
        f.write('B类不确定度: ' + str(b_class_uncertainty) + '\n')
        f.write('不确定度: ' + str(combined_uncertainty) + '\n')
        f.write('相对不确定度: ' + str(100 * combined_uncertainty / mean) + '%\n')
        f.write('\n')
        f.write('测量结果为: \n       x = ' + str(mean) + '±' + str(combined_uncertainty) + '\n')
        f.write('     Δx/x = ' + str(100 * combined_uncertainty / mean) + '%\n')
        

if __name__ == '__main__':
    file_path = 'input.txt' # 修改为你的input.txt文件路径
    data, instrument_error, confidence_coefficient = read_data(file_path)
    good_values, outliers = remove_outliers(data)
    mean = np.mean(good_values)
    std_dev = np.std(good_values, ddof=1)  * t683(good_values)
    save_results(good_values, outliers, mean, std_dev, instrument_error, confidence_coefficient)

测试

在写程序的过程中,计算过程不断的修修补补,东西也比较乱,导致对自己的结果不是很有信心。

拿书上的例题测试了一下,四舍五入后结果正确,但是原始值和书上的值存在一点点不同,应该是这个程序在计算过程中没有进行有效数字处理,微小的偏差积累下来导致有效位数后面的值不准确。再加上 Python 本身浮点计算的误差,导致了最后的结果不太一样。

input.txt

2.125 2.131 2.121 2.127 2.124 2.126
0.004

result.txt

剩余的值: 2.125 2.131 2.121 2.127 2.124 2.126
剔除的值: 
最终的平均数: 2.125666666666667
S = 0.0033266599866331402
σ = 0.0033266599866331402, 3σ = 0.00997997995989942
每个最终值的标准差: 0.0006666666666670373 0.005333333333332746 0.004666666666667041 0.0013333333333327424 0.0016666666666669272 0.00033333333333285253
A类不确定度(标准偏差): 0.001358103252497517
B类不确定度: 0.002309401076758503
不确定度: 0.002679137506321329
相对不确定度: 0.12603751793890522%

测量结果为: //没有减去可定系统误差 -0.008mm
       x = 2.125666666666667 ± 0.002679137506321329
     Δx/x = 0.12603751793890522%

书本上的例题答案

σ = 0.0033
3σ = 0.01
A类不确定度(标准偏差): 0.001
B类不确定度: 0.002
相对不确定度: 0.10%
x = 2.134 ± 0.002 //这一个结果是减去可定系统误差的结果,其实和程序的结果是相同的

后面在实际运用的过程中发现和原有的实验计算工具的不确定度结果差了一个数量级。经过我和同学的反复确认,程序没有问题,原有的工具甚至都没要求输入仪器误差,计算也明显存在问题。

不过,这个程序使用的人次还是偏少,检验的次数也少,需要后面在使用中不断发现问题,改正计算方法,

打包与项目

为了便于更多人使用,在同学的提示下用 Pyinstaller 打包 Python 成可执行文件,最后效果还是不错。

项目也已经上传 Github:

RichardTang2003/physical-experiments-utilities

欢迎star和issues.

For English Reader: A Brief Note on the Development of a Physics Experiment Utility

作者: illlights

2024 © typecho & elise