如何在 python 中将 Quad 实现为列表？

Question

我一直在研究应用一种称为 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, lb, 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,ub,.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,dists,lb,ub)
plt.plot(domain,graph, 'm', label='Conflated Dist.')
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("Conflated PDF")
plt.legend()
plt.show()

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

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,dists,lb,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, lb, ub, 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,dists)/denom
    print('Conflated PDF: ', conflated_pdf)
    return conflated_pdf

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

dist_1_pdf = st.norm.pdf(domain, 2,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, 2.4,1.3)
dist_5_pdf = st.norm.pdf(domain, 2.7,1.5)

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

from matplotlib import pyplot as plt
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("Conflated PDF")
plt.plot(domain, dist_1_pdf, 'r', label='Dist. 1')
plt.plot(domain, dist_2_pdf, 'g', label='Dist. 2')
plt.plot(domain, dist_3_pdf, 'b', label='Dist. 3')
plt.plot(domain, dist_4_pdf, 'y', label='Dist. 4')
plt.plot(domain, dist_5_pdf, 'c', 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_2, dist_3, dist_4, dist_5]
dists=[dist_1_pdf, dist_2_pdf, dist_3_pdf, dist_4_pdf, dist_5_pdf]
graph=conflate_pdf(domain,dists,lb,ub)

plt.plot(domain,graph, 'm', 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,dists):
    p_list=[]
    p_pdf=1
    print('Incoming Array:', p_pdf)
    for dist in dists:
        print('Incoming Distribution Array:', dist.pdf(x))
        p_pdf=p_pdf*dist.pdf(x)
        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

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

def my_product_pdf(x,dists):
    p_list=[]
    p_pdf=1
    print('Incoming Array:', 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_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 my_conflate_pdf(x,dists,lb,ub):
    conflated_pdf_list=[]
    num=my_product_pdf(x,dists)
    print('\n')
    # print('product pdf: ', prod_pdf(x,dists))
    print('product pdf: ', my_product_pdf(x,dists))
    denom = fixed_quad(my_product_pdf, lb, ub, args=(dists,), 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,ub,.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, 2,1)
# dist_2_pdf = st.norm.pdf(domain, 2.5,1.5)
# dist_3_pdf = st.norm.pdf(domain, 2.2,1.6)
# dist_4_pdf = st.norm.pdf(domain, 2.4,1.3)
# dist_5_pdf = st.norm.pdf(domain, 2.7,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)
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)

dist_1_pdf = st.norm.pdf(domain, 2,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, 2.4,1.3)
dist_5_pdf = st.norm.pdf(domain, 2.7,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, dist_1_pdf, 'r', label='Dist. 1')
plt.plot(domain, dist_2_pdf, 'g', label='Dist. 2')
plt.plot(domain, dist_3_pdf, 'b', label='Dist. 3')
plt.plot(domain, dist_4_pdf, 'y', label='Dist. 4')
plt.plot(domain, dist_5_pdf, 'c', label='Dist. 5')

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

# user_graph /= user_graph.sum()

plt.plot(domain, user_graph, 'm', 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, dist_1_pdf, 'r', label='Dist. 1')
plt.plot(domain, dist_2_pdf, 'g', label='Dist. 2')
plt.plot(domain, dist_3_pdf, 'b', label='Dist. 3')
plt.plot(domain, dist_4_pdf, 'y', label='Dist. 4')
plt.plot(domain, dist_5_pdf, 'c', label='Dist. 5')

dists=[dist_1_pdf, dist_2_pdf, dist_3_pdf, dist_4_pdf, dist_5_pdf]
my_graph=my_conflate_pdf(domain,dists,lb,ub)
print('Final Conflated PDF: ', my_graph)

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

# my_graph = inverse_normalise(my_graph)

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

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

这是以下输出：

但是，如果我规范化 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()

输出：

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

Answer 1

第一个大抱怨 - 您没有显示完整的回溯！到现在为止，您应该知道检查回溯以确切了解导致问题的原因是很重要的。

从你的表现来看，我推断问题出在：

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,dists,lb,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 in

denom = quad(product_pdf, lb, 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 需要一个标量值函数。