2016-12-14 62 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爲兩個不同的區域整合。我嘗試了以下的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) 

有沒有更簡單/更清晰的方式去分手了概率密度函數成兩個空間區域,然後整合在每個時間步的每個空間區域內?

回答

1

可以通過使用矢量化操作來消除循環。目前尚不清楚Pd是否爲2D NumPy陣列;它是別的東西(例如,列表列表),它應該被轉換爲一個帶有np.array(...)的2D NumPy數組。之後,你可以這樣做:

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

無需循環的時間索引;切片一次同時發生(具有:代替索引意味着「所有有效索引」)。

同樣,np.trapz可以一次整合所有的時間片:

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

每個P1和P2現在是一個時間序列積分。 axis參數確定沿着哪個軸Pd1進行積分 - 它是第二個軸,即空間。