如何在python中实现Quad到列表？

编程问答 2022-05-14

问题描述

我一直在研究应用一种称为 Conflation 的方法的概念，该方法着眼于将多个分布组合在一起。我在以下链接中问了这个问题：How to implement Conflation for probability distribution in python?，这里用以下代码回答这个问题：

import scipy.stats as st
from scipy.integrate import quad
from scipy import stats
import numpy as np
from matplotlib import pyplot as plt

def prod_pdf(x,dists):
    p_pdf=1
    for dist in dists:
        print('Incoming Array:',dist.pdf(x))
        p_pdf=p_pdf*dist.pdf(x)
        print('final:',p_pdf)
    return print('Product: ',p_pdf),p_pdf
    
def product_pdf(x,dists):
    p_pdf=1
    for dist in dists:
        p_pdf=p_pdf*dist.pdf(x)
    return p_pdf

def conflate_pdf(x,dists,lb,ub):
    print('Input product pdf: ',prod_pdf(x,dists)[1])
    denom = quad(product_pdf,ub,args=(dists,))[0]
    print('Denom: ',denom)
    conflated_pdf=product_pdf(x,dists)/denom
    print('Conflated PDF: ',conflated_pdf)
    return conflated_pdf

lb=-10
ub=10
domain=np.arange(lb,.01)

dist_1 = st.norm.pdf(domain,2,1)
dist_2 = st.norm.pdf(domain,2.5,1.5)
dist_3 = st.norm.pdf(domain,2.2,1.6)
dist_4 = st.norm.pdf(domain,2.4,1.3)
dist_5 = st.norm.pdf(domain,2.7,1.5)

plt.plot(domain,dist_1,'r',label='dist. 1')
plt.plot(domain,dist_2,'g',label='dist. 2')
plt.plot(domain,dist_3,'b',label='dist. 3')
plt.plot(domain,dist_4,'y',label='dist. 4')
plt.plot(domain,dist_5,'c',label='dist. 5')

dists=[stats.norm(2,1),stats.norm(2.5,1.5),stats.norm(2.2,1.6),stats.norm(2.4,1.3),stats.norm(2.7,1.5)]
graph=conflate_pdf(domain,ub)
plt.plot(domain,graph,'m',label='Conflated dist.')
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("Conflated PDF")
plt.legend()
plt.show()

我关注的领域是连续分布部分，它使用分母积分的方程。我知道代码中的错误是由于使用集成时出现 1 大小的数组错误。所以，我认为四边形的输入是使用单个整数，我尝试修改代码以实现与上述代码相同的结果。代码可以在下面找到：

import scipy.stats as st
from scipy.integrate import quad
from scipy import stats
import numpy as np
from matplotlib import pyplot as plt

def product_pdf(x,dists):
    p_list=[]
    p_pdf=1
    print('Incoming Array:',p_pdf)
    list_full_size=np.array(dists).shape
    print('list size: ',list_full_size[0])
    for x in range(list_full_size[1]):
        p_pdf=1
        for y in range(list_full_size[0]):
            p_pdf=p_pdf*dists[y][x]
            # print('Incoming distribution Array:',dist)
            print('Product value: ',p_pdf)
        print('Product PDF:',p_pdf)
        p_list.append(p_pdf)
    print('final Product PDF:',p_pdf)
    print('Product PDF list: ',p_list)
    # return p_pdf
    return p_list
    # return np.array(p_list)

def conflate_pdf(x,ub):
    denom_list=[]
    # time.sleep(500)
    print('\n')
    print('product pdf: ',product_pdf(x,dists))
    prod=product_pdf(x,dists)
    prod_size=np.array(prod).shape[0]
    print('Product Size: ',prod_size)
    # time.sleep(20)
    for i in range(prod_size):
        denom = quad(product_pdf,args=(prod[i],))[0]
        denom_list.append(denom)
        print('Denom: ',denom)
    print('Denominator list: ',denom_list)
    # time.sleep(20)
    # conflated_pdf=prod_pdf(x,dists)/denom
    conflated_pdf=product_pdf(x,.01)

dist_1_pdf = st.norm.pdf(domain,1)
dist_2_pdf = st.norm.pdf(domain,4,2)
dist_3_pdf = st.norm.pdf(domain,7,4)
dist_4_pdf = st.norm.pdf(domain,1.3)
dist_5_pdf = st.norm.pdf(domain,1.5)

# dist_1 = list(st.norm.pdf(domain,1))
# dist_2 = list(st.norm.pdf(domain,1.5))
# dist_3 = list(st.norm.pdf(domain,1.6))
# dist_4 = list(st.norm.pdf(domain,1.3))
# dist_5 = list(st.norm.pdf(domain,1.5))

from matplotlib import pyplot as plt
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("Conflated PDF")
plt.plot(domain,dist_1_pdf,dist_2_pdf,dist_3_pdf,dist_4_pdf,dist_5_pdf,label='dist. 5')

dist_1 = st.norm(2,1)
dist_2 = st.norm(4,2)
dist_3 = st.norm(7,4)
dist_4 = st.norm(2.4,1.3)
dist_5 = st.norm(2.7,1.5)

# dists=[dist_1,dist_5]
dists=[dist_1_pdf,dist_5_pdf]
graph=conflate_pdf(domain,ub)

plt.plot(domain,label='Conflated PDF')
plt.legend()
plt.show()

当我运行代码时，我得到以下输出：

print('list size: ',list_full_size[0])

IndexError: tuple index out of range

我知道错误与列表的大小有关，但是，我该如何修改代码才能产生相同的结果并绘制到第一个代码？

编辑 1：

我再次查看代码并设法发现，如果我将 quad 替换为 fixed_quad，它将能够毫无错误地运行代码。但是，绘图与 user_conflation_pdf 方法不同。这是我的方法的代码，

import scipy.stats as st
import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler,normalizer,normalize,StandardScaler
from scipy.integrate import quad,simps,quad_vec,nquad,cumulative_trapezoid
from scipy.integrate import romberg,trapezoid,simpson,romb
from scipy.integrate import fixed_quad,quadrature,quad_explain
from scipy import stats
import time

def user_prod_pdf(x,p_pdf)
    for dist in dists:
        print('Incoming distribution Array:',dist.pdf(x))
        p_pdf=p_pdf*dist.pdf(x)
        print('Product PDF:',p_list)
    return p_pdf

def user_conflate_pdf(x,user_prod_pdf(x,dists))
    denom = quad(user_prod_pdf,denom)
    conflated_pdf=user_prod_pdf(x,conflated_pdf)
    return conflated_pdf

def my_product_pdf(x,p_pdf)
    list_full_size=np.array(dists).shape
    print('Full list size: ',list_full_size)
    print('list size: ',list_full_size[0])
    for x in range(list_full_size[1]):
        p_pdf=1
        for y in range(list_full_size[0]):
            p_pdf=float(p_pdf)*dists[y][x]
            print('Product value: ',p_list)
    # return p_pdf
    return p_list
    # return np.array(p_list)
    
def my_conflate_pdf(x,ub):
    conflated_pdf_list=[]
    num=my_product_pdf(x,dists)
    print('\n')
    # print('product pdf: ',dists))
    print('product pdf: ',my_product_pdf(x,dists))
    denom = fixed_quad(my_product_pdf,),n=1)[0]
    print('Denom: ',denom)
    # conflated_pdf=prod_pdf(x,dists)/denom
    # conflated_pdf=my_product_pdf(x,dists)/denom
    # conflated_pdf=[i / j for i,j in zip(my_product_pdf(x,dists),denom)]
    for i in np.arange(np.array(num).shape[0]):
        conflated_pdf = num[i] / denom
        conflated_pdf_list.append(conflated_pdf)
    print('Conflated PDF: ',conflated_pdf_list)
    return conflated_pdf_list

lb=-10
ub=10
domain=np.arange(lb,.01)

# dist_1 = st.norm(2,1)
# dist_2 = st.norm(2.5,1.5)
# dist_3 = st.norm(2.2,1.6)
# dist_4 = st.norm(2.4,1.3)
# dist_5 = st.norm(2.7,1.5)

# dist_1_pdf = st.norm.pdf(domain,1)
# dist_2_pdf = st.norm.pdf(domain,1.5)
# dist_3_pdf = st.norm.pdf(domain,1.6)
# dist_4_pdf = st.norm.pdf(domain,1.3)
# dist_5_pdf = st.norm.pdf(domain,1.5)

# dist_1_pdf /= dist_1_pdf.sum()
# dist_2_pdf /= dist_2_pdf.sum()
# dist_3_pdf /= dist_3_pdf.sum()
# dist_4_pdf /= dist_4_pdf.sum()
# dist_5_pdf /= dist_5_pdf.sum()

dist_1 = st.norm(2,1.5)

dist_1_pdf = st.norm.pdf(domain,1.5)

# dist_1_pdf /= dist_1_pdf.sum()
# dist_2_pdf /= dist_2_pdf.sum()
# dist_3_pdf /= dist_3_pdf.sum()
# dist_4_pdf /= dist_4_pdf.sum()
# dist_5_pdf /= dist_5_pdf.sum()

# User:
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("User Conflated PDF")
plt.plot(domain,label='dist. 5')

dists=[dist_1,dist_5]
user_graph=user_conflate_pdf(domain,ub)
print('Final Conflated PDF: ',user_graph)

# user_graph /= user_graph.sum()

plt.plot(domain,user_graph,label='Conflated PDF')
plt.legend()
plt.show()

# My Code:
# from matplotlib import pyplot as plt
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("My Conflated PDF Code")
plt.plot(domain,label='dist. 5')

dists=[dist_1_pdf,dist_5_pdf]
my_graph=my_conflate_pdf(domain,my_graph)

# my_graph /= np.array(my_graph).sum()

# my_graph = inverse_normalise(my_graph)

plt.plot(domain,my_graph,label='Conflated PDF')
plt.legend()
plt.show()

# Conflated PDF:
print('User Conflated PDF: ',user_graph)
print('My Conflated PDF: ',my_graph)

输出如下：

enter image description here

enter image description here

但是，如果我规范化 pdf 列表，输出图是相同的。以下是需要更改以规范列表的内容，

user_graph /= user_graph.sum()
dist_1_pdf /= dist_1_pdf.sum()
dist_2_pdf /= dist_2_pdf.sum()
dist_3_pdf /= dist_3_pdf.sum()
dist_4_pdf /= dist_4_pdf.sum()
dist_5_pdf /= dist_5_pdf.sum()

输出：

enter image description here

enter image description here

如何修改我的代码/函数才能获得此图？

enter image description here

解决方法

第一个大抱怨 - 您没有显示完整的追溯！现在您应该知道检查回溯以准确找出导致问题的原因是必不可少的。

从你所展示的一点点来看，我推断问题出在：

def product_pdf(x,dists):
    p_list=[]
    p_pdf=1
    print('Incoming Array:',p_pdf)
    list_full_size=np.array(dists).shape
    print('list size: ',list_full_size[0])
    ...

你为什么不突出那个地方？阅读代码需要花费不必要的时间才能找到打印行。

该函数在 conflate_pdf(x,dists,lb,ub) 中被多次调用。至少一次打印。那个打印成功了吗？您的代码有一堆打印件，但您没有显示任何打印件:(

来自

dists=[dist_1_pdf,dist_2_pdf,dist_3_pdf,dist_4_pdf,dist_5_pdf]
graph=conflate_pdf(domain,ub)

dists 是 5 个对象的列表。因此，我希望 np.array(dists) 至少是一个一维数组，如果那些 pdf 对象是具有相同形状的数组，则可能更多。这意味着 list_full_size[0] 应该可以工作，并且匹配 len(dists)（即不需要将列表转换为数组。

array 转换排序有助于

for x in range(list_full_size[1]):

虽然也可以

for x in range(len(dists[0])):

因为 x 和 y 循环从 dists 而不是 np.array(dists) 工作。

无论如何，看起来 product_pdf(x,dists) 调用应该可以工作。

但是，您也在

中使用了 product_pdf 
denom = quad(product_pdf,ub,args=(prod[i],))[0]

dists 参数现在是 prod[i]，而不是原来的 dists。

什么是prod？

prod=product_pdf(x,dists)

等等，product_pdf 不是应该为 quad 生成一个的数字吗？你用 prod[i] 做什么，并将它作为 dist 参数传递给 product_pdf。

product_pdf 返回 plist，它看起来像一个数字列表。所以 prod[i] 是一个数字，它的形状，作为一个数组是 ()。因此错误。现在我们又回到了我上次帮助您解决的错误 - quad 需要一个标量值函数。

arraysarraysarrayslistlistnumpypythonscipyscipy