2016-12-14 15 views
1

私はいくつかの確率密度関数を持っている:分割統合確率密度

T = 10000 
tmin = 0 
tmax = 10**20 
t = np.linspace(tmin, tmax, T) 
time = np.asarray(t)     #this line may be redundant 
for j in range(T): 
    timedep_PD[j]= probdensity_func(x,time[j],initial_state) 

私はxの2つの異なる領域の上にそれを統合したいです。私は2つの空間領域にtimedep_PD配列を分割するために、次の試み、その後、統合を進め:

step = abs(xmin - xmax)/T 
l1 = int(np.floor((abs(ab - xmin)* T)/abs(xmin - xmax))) 
l2 = int(np.floor((abs(bd - ab)* T)/abs(xmin - xmax))) 

#For spatial region 1 
R1 = np.empty([l1]) 
R1 = x[:l1] 
for i in range(T): 
    Pd1[i] = Pd[i][:l1] 

#For spatial region 2 
Pd2 = np.empty([T,l2]) 
R2 = np.empty([l2]) 
R2 = x[l1:l1+l2] 
for i in range(T): 
    Pd2[i] = Pd[i][l1:l1+l2] 

#Integrating over each spatial region 
for i in range(T): 
    P[0][i] = np.trapz(Pd1[i],R1) 
    P[1][i] = np.trapz(Pd2[i],R2) 

は、2つの空間領域に確率密度関数を分割して、統合について移動する簡単/より明確な方法はあります各時間領域の各空間領域内で?

答えて

1

ループを削除するには、ベクトル化された操作を使用します。 Pdが2D NumPy配列であるかどうかは不明です。それは何か他のもの(リストのリストなど)で、np.array(...)で2D NumPy配列に変換する必要があります。その後、これを行うことができます:

Pd1 = Pd[:, :l1] 
Pd2 = Pd[:, l1:l1+l2] 

時間インデックスをループする必要はありません。一度にスライシングが行われます(インデックスの代わりに:は「すべての有効なインデックス」を意味します)。

P1 = np.trapz(Pd1, R1, axis=1) 
P2 = np.trapz(Pd2, R2, axis=1) 

それぞれP1及びP2は、現在の積分の時系列である:

同様に、np.trapzは、一度にすべてのタイムスライスを統合することができます。 axisパラメータは、どの軸に沿ってPd1が統合されるかを決定します。これは、第2の軸、つまり空間です。

関連する問題