如何在 pandas 数据框中进行不确定性传播（地球化学数据缩减）

Question

我尝试根据我的标准校准曲线计算其化学浓度。

我有四天的数据，每一个都有自己的标准曲线。我已经对标准进行了约克回归，并获得了以下截距 (a)、斜率 (b) 和相关误差矩阵：

york = {'date':['Jun27', 'Jun28', 'Jun29', 'Jun30'],
        'a':[1.2013, 1.0057, 1.1462, 0.3874],
        'b':[44138,41246,43311,49830],
        'siga':[0.2795,0.2791,0.2189,0.3641],
        'sigb':[531.7,873.7,727.26,1251.7]}

yk = pd.DataFrame(york) 
yk.set_index('date', inplace = True)

使得[Ti] = 49Ti/30Si16O * b + a

我也有数据；通常我使用

将它作为数据框读入

df30 = pd.read_clipboard()

因为它是一个很大的块，这样我就可以从复制的电子表格中获取列名。但举个例子，这里有几个数据点：

Jun30 = {'File':['LB13-LP41-10-ZR.asc', 'LB13-LP41-19-ZR.asc', 'LB13-MB50-1-ZR.asc', 'LB13-MB50-18ZR.asc'],
         '49Ti/30Si16O':[0.000405567, 0.000272094, 0.000320981, 0.000153742],
         '1 se err':[2.61586E-06, 7.65216E-07, 1.32338E-06, 1.53561E-06]}
        df30 = pd.DataFrame(Jun30)
        df30.set_index('File', inplace = True)

我想对分析误差和标准校准误差进行蒙特卡罗不确定性传播，这样

[Ti]+/- [Tierr] = (49Ti/30Si16O+/-1 se err) * (b +/- sigb) + (a +/- siga)

在数据框中执行此操作的最简单/最经济的方法是什么？理想情况下，我想在数据框中添加两列：“[Ti]”和“Ti err”，但我不知道如何遍历每一行并引用正确的值。

通常在 MATLAB 中使用数组执行此操作，我会执行以下操作：

    RTi = [data for Ti ratio]
    RTierr = [associated errors]
    %etc...
    N = numel(RTi)
    Ti = zeros(N,1);
    Tierr = zeros(N,1);

    for i = 1:N
        j = zeros(1e5,1);
        k = zeros(1e5,1);

        for n = 1:1e5
        a(n) = normrnd(intercept,sigintercept);
        b(n) = normrnd(slope,sigslope);
        k(n) = normrnd(RTi,RTierr);
        j(n) = k(n).*b(n)+a(n)
    end
    Ti(i) = mean(j);
    Tierr(i) = std(j);
end

但这有点笨拙，我很确定学习如何在数据帧中使用 python 做到这一点会更容易，希望更快。

Answer 1

给定以下形式的数据：

Jun27 = {'File':['LB13-LP41-10-ZR.asc', 'LB13-LP41-19-ZR.asc', 'LB13-MB50-1-ZR.asc', 'LB13-MB50-18ZR.asc'],
         '49Ti/30Si16O':[0.000405567, 0.000272094, 0.000320981, 0.000153742],
         '1 se err':[2.61586E-06, 7.65216E-07, 1.32338E-06, 1.53561E-06], 'date': 'Jun27'}
Jun28 = {'File':['LB13-LP41-10-ZR.asc', 'LB13-LP41-19-ZR.asc', 'LB13-MB50-1-ZR.asc', 'LB13-MB50-18ZR.asc'],
         '49Ti/30Si16O':[0.000405567, 0.000272094, 0.000320981, 0.000153742],
         '1 se err':[2.61586E-06, 7.65216E-07, 1.32338E-06, 1.53561E-06], 'date': 'Jun28'}
Jun29 = {'File':['LB13-LP41-10-ZR.asc', 'LB13-LP41-19-ZR.asc', 'LB13-MB50-1-ZR.asc', 'LB13-MB50-18ZR.asc'],
         '49Ti/30Si16O':[0.000405567, 0.000272094, 0.000320981, 0.000153742],
         '1 se err':[2.61586E-06, 7.65216E-07, 1.32338E-06, 1.53561E-06], 'date': 'Jun29'}
Jun30 = {'File':['LB13-LP41-10-ZR.asc', 'LB13-LP41-19-ZR.asc', 'LB13-MB50-1-ZR.asc', 'LB13-MB50-18ZR.asc'],
         '49Ti/30Si16O':[0.000405567, 0.000272094, 0.000320981, 0.000153742],
         '1 se err':[2.61586E-06, 7.65216E-07, 1.32338E-06, 1.53561E-06], 'date': 'Jun30'}

• 只提供了一天的数据，因此在本示例中，所有天都使用了这一天

连接所有数据：

data = pd.concat([pd.DataFrame(Jun27), pd.DataFrame(Jun28), pd.DataFrame(Jun29), pd.DataFrame(Jun30)])

• 注意添加的日期列

合并yk（来自示例）和data：

df = pd.merge(yk, data)

• 使用同一DataFrame 中的所有数据执行计算更容易
• 合并在date 列上，同时包含在DataFrames 中

创建计算：

df['Ti'] = df['49Ti/30Si16O'] * df.b + df.a

df['49Ti/30Si16O_error_min'] = df['49Ti/30Si16O'] - df['1 se err']
df['49Ti/30Si16O_error_max'] = df['49Ti/30Si16O'] + df['1 se err']
df['b_error_min'] = df.b - df.sigb
df['b_error_max'] = df.b + df.sigb
df['a_error_min'] = df.a - df.siga
df['a_error_max'] = df.a + df.siga
df['Ti_min'] = df['49Ti/30Si16O_error_min'] * df['b_error_min'] + df['a_error_min']
df['Ti_max'] = df['49Ti/30Si16O_error_max'] * df['b_error_max'] + df['a_error_max']