如何在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误差）*（b+/-sigb）+（a+/-siga）

在数据帧中，最简单/最经济的方法是什么？理想情况下，我想在dataframe中添加两列：“[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)])

• 注意添加了date列

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

df = pd.merge(yk, data)

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

创建计算：

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']