【问题标题】:Quadratic n term equation using multiindex使用多索引的二次 n 项方程
【发布时间】:2019-04-25 13:33:12
【问题描述】:

我有两个 DF,我想用它们来计算以下内容:

w(ti,ti)*a(ti)^2 + w(tj,tj)*b(sj,tj)^2 + 2*w(si,tj)*a(ti)*b(tj)

上面使用了两个术语(a,b)。 w 是权重 df,其中 i 和 j 是属于 a 和 b 的 Tn 索引的索引和列空间。

设置 - 编辑动态 W

import pandas as pd
import numpy as np

I = ['i'+ str(i) for i in range(4)]
Q = ['q' + str(i) for i in range(5)]
T = ['t' + str(i) for i in range(3)]
n = 100

df1 = pd.DataFrame({'I': [I[np.random.randint(len(I))] for i in range(n)],
                    'Q': [Q[np.random.randint(len(Q))] for i in range(n)],
                    'Tn': [T[np.random.randint(len(T))] for i in range(n)],
                    'V': np.random.rand(n)}).groupby(['I','Q','Tn']).sum()

df1.head(5)
I  Q  Tn  V        
i0 q0 t0  1.626799
      t2  1.725374
   q1 t0  2.155340
      t1  0.479741
      t2  1.039178

w = np.random.randn(len(T),len(T))
w = (w*w.T)/2
np.fill_diagonal(w,1)
W = pd.DataFrame(w, columns = T, index = T)

W
          t0        t1        t2
t0  1.000000  0.029174 -0.045754
t1  0.029174  1.000000  0.233330
t2 -0.045754  0.233330  1.000000

实际上,我想使用 df1 中的索引 Tn 来为每个 I 和 Q 使用上述等式。

上例中df1.loc['i0','q0'] 的最终结果应该是:

  W(t0,t0) * V(t0)^2 
+ W(t2,t2) * V(t2)^2
+ 2 * W(t0,t2) * V(t0) * V(t2) 
=     
  1.0 * 1.626799**2 
+ 1.0 * 1.725374**2 
+ (-0.045754) * 1.626799 * 1.725374

上例中df1.loc['i0','q1'] 的最终结果应该是:

  W(t0,t0) * V(t0)^2 
+ W(t1,t1) * V(t1)^2
+ W(t2,t2) * V(t2)^2
+ 2 * W(t0,t1) * V(t0) * V(t1)
+ 2 * W(t0,t2) * V(t0) * V(t2)
+ 2 * W(t2,t1) * V(t1) * V(t2)
=     
  1.0 * 2.155340**2 
+ 1.0 * 0.479741**2
+ 1.0 * 1.039178**2
+ 0.029174 * 2.155340 * 0.479741 * 1
+ (-0.045754) * 2.155340 * 1.039178 * 1
+ 0.233330 * 0.479741 * 1.039178 * 1

此模式将根据每个 Q 中的 tn 项的数量重复,因此它应该足够健壮,可以根据需要处理尽可能多的 Tn 项(在示例中我使用 3,但它可能多达 100 个或更多) .

然后应将每个结果保存在带有Index = [I, Q] 的新 DF 中 当n 值增加时,该解决方案也不应该比 excel 慢。

提前致谢

【问题讨论】:

  • 你的方程意味着值'w'对于所有三个项都是相同的,但它们不是。也许您应该重命名它们并描述它们与 df1 indices 的关系或派生方式。让您的读者更轻松。
  • df1.loc['i0','q0' 有三个Tn。它是如何工作的?
  • W 不应该是对称的吗?如果不是,对于您给出的示例,我怎么知道在W.loc['t3','t4']W.loc['t4','t3'] 之间使用哪个因素,因为您使用的是第一个因素,但为什么呢?
  • 我已将问题更改为与 cmets 对应

标签: python numpy dataframe multi-index quadratic


【解决方案1】:

一种方法可能是首先reindex您的数据框df1与列表IQTnpd.MultiIndex.from_product的所有可能组合,填充“V”列中的缺失值为 0。然后该列具有 len(I)*len(Q)*len(T) 元素。然后您可以通过reshapevalues 获取与IQ 上的一个组合相关的每一行,例如:

ar = (df1.reindex(pd.MultiIndex.from_product([I,Q,T], names=['I','Q','Tn']),fill_value=0)
         .values.reshape(-1,len(T)))

要查看我的输入 df1ar 之间的关系,这里有一些相关的行

print (df1.head(6))
                 V
I  Q  Tn          
i0 q0 t1  1.123666
   q1 t0  0.538610
      t1  2.943206
   q2 t0  0.570990
      t1  0.617524
      t2  1.413926
print (ar[:3])
[[0.         1.1236656  0.        ]
 [0.53861027 2.94320574 0.        ]
 [0.57099049 0.61752408 1.4139263 ]]

现在,要与W 的元素相乘,一种方法是创建ar 与自身的外积,但逐行获取,每行一个len(T)*len(T) 矩阵。例如,对于第二行:

[0.53861027 2.94320574 0.        ]

变成

[[0.29010102, 1.58524083, 0.        ], #0.29010102 = 0.53861027**2, 1.58524083 = 0.53861027*2.94320574 ...
 [1.58524083, 8.66246003, 0.        ],
 [0.        , 0.        , 0.        ]]

可能有几种方法,例如ar[:,:,None]*ar[:,None,:]np.einsum 右下标:np.einsum('ij,ik->ijk',ar,ar)。两者都给出相同的结果。

下一步可以使用tensordot 并指定正确的axes。因此,使用 arW 作为输入,您可以:

print (np.tensordot(np.einsum('ij,ik->ijk',ar,ar),W.values,axes=([1,2],[0,1])))
array([ 1.26262437, 15.29352438, 15.94605435, ...

在这里检查第二个值,1*0.29010102 + 1*8.66246003 + 2.*2*1.58524083 == 15.29352438(其中 1 是 W(t0,t0)W(t1,t1),2 是 W(t0,t1)

最后,按预期创建数据框,再次使用pd.MultiIndex.from_product

new_df = pd.DataFrame({'col1': np.tensordot(np.einsum('ij,ik->ijk',ar,ar),
                                            W.values,axes=([1,2],[0,1]))},
                      index=pd.MultiIndex.from_product([I,Q], names=['I','Q']))

print (new_df.head(3))
            col1
I  Q            
i0 q0   1.262624
   q1  15.293524
   q2  15.946054
...

注意:如果你确定T的每个元素至少在df1的最后一级出现一次,则ar可以使用unstack如@ 987654359@。但我建议使用上面的reindex 方法来防止任何错误

【讨论】:

  • 这似乎有效。但是,我在我的问题中发现了一个边缘情况,这会使这个答案不正确。否则你教会了我一些新东西!谢谢
  • @RealRageDontQuit 你所谓的边缘情况实际上是一个不同的问题。您可以通过添加索引级别来更改数据框结构,通过与另一个矩阵 s 相乘来更改公式,并对这个新索引进行求和(至少我所理解的)。我认为我的回答可以很容易地适应这个问题,但如果你想要一个通用的方法,它会更复杂。
  • 你说这个问题解决了最初的问题是正确的。我已更改问题以反映最初的问题,将您的答案勾选为已接受(谢谢),并且还将使用这种极端情况创建一个新问题。
猜你喜欢
  • 2018-06-20
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2020-11-28
  • 1970-01-01
  • 1970-01-01
  • 2019-04-13
  • 1970-01-01
相关资源
最近更新 更多