问题描述
我一直在考虑生成一个乘法函数，用于一种称为 Conflation 的方法。该方法可以在以下文章 (An Optimal Method for Consolidating Data from Different Experiments) 中找到。可以在下面找到融合方程：
[a*b for a,b in zip(lista,listb)]
list(map(operator.mul,lista,listb))
np.multiply(lista,listb)
ab = [lista[i]*listb[i] for i in range(len(lista))]
lista = [1,2,3,4]
listb = [2,4,5]
ab = [] #Create empty list
for i in range(0,len(lista)):
ab.append(lista[i]*listb[i]) #Adds each element to the list
但是查看 2 个以上的列表，我不断收到有关大小为 1 数组的错误消息，或者代码查看每个分布中的第一个变量，但是，对于循环的其余部分，它继续打印相同的值，它不会转到列表中的下一个值，并且 Conflated 分布是单个变量。请参阅以下代码以及部分输出和错误消息：
from scipy.integrate import quad
from scipy import stats
import numpy as np
def prod_pdf(x,dists):
p_pdf=1
print('Incoming Array:',p_pdf)
for c,dist in enumerate(dists):
p_pdf=p_pdf*dist[c]
print('final:',p_pdf)
return p_pdf
def conflate_pdf(x,dists,lb,ub):
print('Input product pdf: ',prod_pdf(x,dists))
denom = quad(prod_pdf,ub,args=(dists,))[0]
# denom = simps(prod_pdf)
# denom = nquad(func=(prod_pdf),ranges=([lb,ub]),))[0]
print('Denom: ',denom)
conflated_pdf=prod_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,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)
from matplotlib import pyplot as plt
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("Conflated PDF")
plt.legend()
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=[dist_1,dist_5]
print('distribution list: \n',dists)
graph=conflate_pdf(domain,ub)
plt.plot(domain,graph,'m',label='Conflated dist.')
plt.show()
输出的一部分：
Incoming Array: 1
final: 2.1463837356630605e-32
final: 5.0231307782193034e-48
final: 3.266239495519432e-61
final: 2.187514996217005e-81
final: 1.979657878680375e-97
Incoming Array: 1
final: 2.1463837356630605e-32
final: 5.0231307782193034e-48
final: 3.266239495519432e-61
final: 2.187514996217005e-81
final: 1.979657878680375e-97
Denom: 3.95931575736075e-96
Incoming Array: 1
final: 2.1463837356630605e-32
final: 5.0231307782193034e-48
final: 3.266239495519432e-61
final: 2.187514996217005e-81
final: 1.979657878680375e-97
Conflated PDF: 0.049999999999999996
第二个代码：
import winsound
from functools import reduce
from itertools import chain
import scipy.stats as st
from glob import glob
from collections import defaultdict,Counter
from sklearn.neighbors import KDTree
import pywt
import peakutils
import scipy
import os
from scipy import signal
from scipy.fftpack import fft,fftfreq,rfft,rfftfreq,dst,idst,dct,idct
from scipy.signal import find_peaks,find_peaks_cwt,argrelextrema,welch,lfilter,butter,savgol_filter,medfilt,freqz,filtfilt
from pylab import *
import glob
import sys
import re
from numpy import NaN,Inf,arange,isscalar,asarray,array
from scipy.stats import skew,kurtosis,median_absolute_deviation
import warnings
import numpy as np
import pandas as pd
import scipy.stats as st
import matplotlib.pyplot as plt
from scipy.stats import pearsonr,kendalltau,spearmanr,ppcc_max
import matplotlib.mlab as mlab
from statsmodels.graphics.tsaplots import plot_acf
from tsfresh.feature_extraction.feature_calculators import mean_abs_change as mac
from tsfresh.feature_extraction.feature_calculators import mean_change as mc
from tsfresh.feature_extraction.feature_calculators import mean_second_derivative_central as msdc
from pyAudioAnalysis.ShortTermFeatures import energy as stEnergy
import pymannkendall as mk_test
from sklearn.preprocessing import MinMaxScaler,normalizer,normalize,StandardScaler
from scipy.integrate import quad,simps,quad_vec,nquad
def prod_pdf(x,dists):
i=0
# p_pdf=np.ones(np.array(dists)[0].shape)
dist_size = np.array(dists).shape
print('Incoming Array:',dists)
print('Incoming Array Size:',dist_size[1])
print('Full Incoming Array Size:',dist_size)
print('Number of Incoming Array Size:',dist_size[0])
# print('Incoming Product Array:',p_pdf)
# print('Incoming Product Array Size:',np.array(p_pdf).shape)
if dist_size[0]==2:
p_pdf=dists[0]*dists[1]
print('final:',p_pdf)
results=dists[0]*dists[1]
i+=1
elif dist_size[0]>2:
results=dists[0]*dists[1]
for i in range(2,dist_size[0]):
p_pdf=results*dists[i]
print('final:',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.legend()
plt.plot(domain,dist_5]
# print('distribution list: \n',label='Conflated dist.')
plt.show()
错误信息：
in line 79,in conflate_pdf:
denom = quad(prod_pdf,))[0]
File "D:\Anaconda\lib\site-packages\scipy\integrate\quadpack.py",line 351,in quad
retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,File "D:\Anaconda\lib\site-packages\scipy\integrate\quadpack.py",line 463,in _quad
return _quadpack._qagse(func,limit)
TypeError: only size-1 arrays can be converted to Python scalars
在我看来，第一个代码是我想要的方法，第二个代码中的错误是因为集成部分需要一个标量数。如何修复这两个代码以获得以下输出？
代码：
from scipy.integrate import quad
from scipy import stats
import numpy as np
def prod_pdf(x,p_pdf)
for dist in dists:
p_pdf=p_pdf*dist.pdf(x)
print('final:',.01)
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)]
from matplotlib import pyplot as plt
plt.plot(domain,label='dist. 5')
graph=conflate_pdf(domain,ub)
plt.plot(domain,label='Conflated dist.')
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("Conflated PDF")
plt.legend()
plt.show()
这是所需输出的一小部分：
Incoming Array: 1
final: 0.15352177537004433
final: 0.034348669264845304
final: 0.006519131844904635
final: 0.0015040030811035296
final: 0.0003607258742065213
Incoming Array: 1
final: 0.042345986284209325
final: 0.006294747321619583
final: 0.0007651214249593444
final: 9.805307029794648e-05
final: 1.668121592516301e-05
Denom: 0.0029066671327537714
Incoming Array: 1
final: [2.14638374e-32 2.41991991e-32 2.72804284e-32 ... 6.41980576e-15
5.92770938e-15 5.47278628e-15]
final: [4.75178372e-48 5.66328097e-48 6.74864868e-48 ... 7.03075979e-21
6.27970218e-21 5.60806584e-21]
final: [2.80912097e-61 3.51131870e-61 4.38823989e-61 ... 1.32670185e-26
1.14952951e-26 9.95834610e-27]
final: [1.51005552e-81 2.03116529e-81 2.73144352e-81 ... 1.76466623e-34
1.46198598e-34 1.21092834e-34]
final: [1.09076800e-97 1.55234627e-97 2.20861552e-97 ... 3.72095218e-40
2.98464396e-40 2.39335035e-40]
Conflated PDF: [3.75264162e-95 5.34063998e-95 7.59844666e-95 ... 1.28014389e-37
1.02682689e-37 8.23400219e-38]
想要的剧情：
编辑 1：
我设法根据@MaxPierini 更新了更新代码，但是，我无法获得所需的混合分布图。查看以下代码和输出：
代码：
import winsound
from functools import reduce
from itertools import chain
import scipy.stats as st
from glob import glob
from collections import defaultdict,dists):
p_pdf=np.ones(np.array(dists)[0].shape)
# p_pdf=1
print('Incoming Array:',np.array(dists)[1].shape)
print('Incoming Product Array:',p_pdf)
print('Incoming Product Array Size:',np.array(p_pdf).shape)
for c,p_pdf)
return p_pdf
# def conflate_pdf(x,ub):
# print('Input product pdf: ',dists))
# denom = quad(prod_pdf,))[0]
# # denom = simps(prod_pdf)
# # denom = nquad(func=(prod_pdf),))[0]
# print('Denom: ',denom)
# conflated_pdf=prod_pdf(x,dists)/denom
# print('Conflated PDF: ',conflated_pdf)
# return conflated_pdf
# use computed PDFs and matrix
def conflate_pdf(x,dists):
# numerator (product)
# num = np.array(dists).prod(axis=0)
num = prod_pdf(x,dists)
print('Input product pdf: ',num)
# conflation = prod_pdf(x,dists)
# normalize (integral)
conflated_pdf = num / num.sum()
print('Conflated PDF: ',1)
# dist_1 /= dist_1.sum()
dist_2 = st.norm.pdf(domain,1.5)
# dist_2 /= dist_2.sum()
dist_3 = st.norm.pdf(domain,1.6)
# dist_3 /= dist_3.sum()
dist_4 = st.norm.pdf(domain,1.3)
# dist_4 /= dist_4.sum()
dist_5 = st.norm.pdf(domain,1.5)
# dist_5 /= dist_5.sum()
# dist_1 = list(st.norm.pdf(domain,dists)
# graph=conflate_pdf(domain,ub)
graph=conflate_pdf(domain,dists)
plt.plot(domain,label='Conflated dist.')
plt.show()
输出：
Incoming Array: [array([2.14638374e-32,2.41991991e-32,2.72804284e-32,...,6.41980576e-15,5.92770938e-15,5.47278628e-15]),array([2.21385563e-16,2.34027620e-16,2.47380598e-16,1.09516706e-06,1.05938091e-06,1.02471859e-06]),array([5.91171893e-14,6.20014921e-14,6.50239789e-14,1.88699641e-06,1.83054781e-06,1.77571847e-06]),array([5.37554463e-21,5.78462242e-21,6.22446263e-21,1.33011515e-08,1.27181248e-08,1.21599343e-08]),array([7.22336360e-17,7.64263883e-17,8.08589121e-17,2.10858694e-06,2.04149972e-06,1.97645911e-06])]
Incoming Array Size: (2000,)
Incoming Product Array: [1. 1. 1. ... 1. 1. 1.]
Incoming Product Array Size: (2000,)
final: [2.14638374e-32 2.14638374e-32 2.14638374e-32 ... 2.14638374e-32
2.14638374e-32 2.14638374e-32]
final: [5.02313078e-48 5.02313078e-48 5.02313078e-48 ... 5.02313078e-48
5.02313078e-48 5.02313078e-48]
final: [3.2662395e-61 3.2662395e-61 3.2662395e-61 ... 3.2662395e-61 3.2662395e-61
3.2662395e-61]
final: [2.187515e-81 2.187515e-81 2.187515e-81 ... 2.187515e-81 2.187515e-81
2.187515e-81]
final: [1.97965788e-97 1.97965788e-97 1.97965788e-97 ... 1.97965788e-97
1.97965788e-97 1.97965788e-97]
Input product pdf: [1.97965788e-97 1.97965788e-97 1.97965788e-97 ... 1.97965788e-97
1.97965788e-97 1.97965788e-97]
Conflated PDF: [0.0005 0.0005 0.0005 ... 0.0005 0.0005 0.0005]
剧情：
编辑 2：
我实现了答案中的代码（由 @MaxPierini 提供），它似乎有效，而且，我设法用 quad 解决了问题，如果我更改了quad 到 fixed_quad 并规范化 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,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,ub):
print('Input product pdf: ',user_prod_pdf(x,dists))
denom = quad(user_prod_pdf,))[0]
print('Denom: ',denom)
conflated_pdf=user_prod_pdf(x,dists)/denom
print('Conflated PDF: ',conflated_pdf)
return conflated_pdf
def user_conflate_pdf_2(pdfs):
"""
Compute conflation of given pdfs.
[ARGS]
- pdfs: PDFs numpy array of shape (n,x)
where n is the number of PDFs
and x is the variable space.
[RETURN]
A 1d-array of normalized conflated PDF.
"""
# conflate
conflation = np.array(pdfs).prod(axis=0)
# normalize
conflation /= conflation.sum()
return conflation
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_pdf)
print('Product PDF:',p_list)
# return p_pdf
return p_list
# return np.array(p_list)
def my_conflate_pdf(x,ub):
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)]
print('Conflated PDF: ',conflated_pdf)
return conflated_pdf
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)
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,1)
dist_2_pdf = st.norm.pdf(domain,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_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,dist_2_pdf,dist_3_pdf,dist_4_pdf,dist_5_pdf,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()
# User 2:
plt.xlabel("domain")
plt.ylabel("pdf")
plt.title("User Conflated PDF 2")
plt.plot(domain,label='dist. 5')
dists=[dist_1_pdf,dist_5_pdf]
user_graph=user_conflate_pdf_2(dists)
print('Final User Conflated PDF 2 : ',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_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: ',np.array(my_graph))
输出如下：
我的问题在这里，我知道我需要规范化 PDF 列表。但是，假设我没有对 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()
我没有标准化的混淆代码图：
解决方法
在第二个 prod_pdf 中，您使用的是计算 PDF，而在第一个中您使用的是定义的分布。因此，在第二个 prod_pdf 中，您已经拥有了 PDF。因此，在 for 循环中，您只需执行 p_pdf = p_pdf * pdf
从您链接的论文中，我们知道“对于离散输入分布，合并的类似定义是概率质量函数的归一化乘积”。因此，您不仅需要获取 PDF 的产品，还需要对其进行规范化。因此，重写离散分布的方程，我们得到
其中 F 是我们需要合并的分布数量，N 是离散变量 x 的长度。
import numpy as np
import scipy.stats as sps
import matplotlib.pyplot as plt
# use defined distributions
def prod_pdf_1(x,dists):
num = 1
for dist in dists:
num *= dist.pdf(x)
den = 0
for i in x:
_mul = 1
for dist in dists:
_mul *= dist.pdf(i)
den += _mul
return num / den
# use computed PDFs
def prod_pdf_2(pdfs):
num = 1
for pdf in pdfs:
num *= pdf
den = 0
for i in range(len(num)):
_mul = 1
for pdf in pdfs:
_mul *= pdf[i]
den += _mul
return num / den
第一个定义使用分布，第二个定义使用 PDF。
现在让我们定义分布和 PDF。
# define x variable
x = np.linspace(-2.5,7.5,100)
# define distributions
dists = [
sps.norm(2.0,1.0),sps.norm(2.5,1.5),sps.norm(2.2,1.6),sps.norm(2.4,1.3),sps.norm(2.7,]
# compute PDFs
pdfs = [
sps.norm.pdf(x,2.0,sps.norm.pdf(x,2.5,2.2,2.4,2.7,]
我们现在可以计算合并和绘图。请注意，我们不需要对混合分布进行归一化，因为我们已经完成了，但我们需要在绘图之前对单个分布进行归一化（总和为 1）。
# first method
p_pdf_1 = prod_pdf_1(x,dists)
# second method
p_pdf_2 = prod_pdf_2(pdfs)
# compare
for pdf in pdfs:
# normalize PDF to sum to 1
pdf /= pdf.sum()
plt.plot(x,pdf)
plt.plot(x,p_pdf_1,label='prod 1',lw=5,color='C1')
plt.plot(x,p_pdf_2,label='prod 2',ls='--',color='k')
plt.legend();
也许这不是最优雅的解决方案，您最好使用矩阵，尤其是当您需要计算大量分布的合并时。
更新
一个更优雅的解决方案是简单地做
# use computed PDFs and matrix
def conflate(pdfs):
# numerator (product)
conflation = pdfs.prod(axis=0)
# normalize (divide by the integral)
conflation /= conflation.sum()
return conflation
并将PDF定义为二维数组（FxN矩阵，F分布PDF和N离散变量x的长度）>
# compute PDFs
pdfs = np.array([
sps.norm.pdf(x,])
我们这样做
conflation = conflate(pdfs)
# compare
for i,pdf in enumerate(pdfs):
# normalize PDF to sum to 1
pdf /= pdf.sum()
plt.plot(x,pdf,label=f'dist {i+1}')
plt.plot(x,conflation,label='conflation',color='k')
plt.legend();
更新 2
完整代码
import numpy as np
import scipy.stats as sps
from matplotlib import pyplot as plt
def conflate(pdfs):
"""
Compute conflation of given pdfs.
[ARGS]
- pdfs: PDFs numpy array of shape (n,x)
where n is the number of PDFs
and x is the variable space.
[RETURN]
A 1d-array of normalized conflated PDF.
"""
# conflate
conflation = pdfs.prod(axis=0)
# normalize
conflation /= conflation.sum()
return conflation
# define x limits and size
lb = -10
ub = 10
size = 1000
# x linear space
x = np.linspace(lb,ub,size)
# define PDFs in x:
# these are Probability Density Functions
# evaluated in x defined linear space
pdf_1 = sps.norm.pdf(x,2,1)
pdf_2 = sps.norm.pdf(x,1.5)
pdf_3 = sps.norm.pdf(x,1.6)
pdf_4 = sps.norm.pdf(x,1.3)
pdf_5 = sps.norm.pdf(x,1.5)
# PDFs in (n,x) array
pdfs = np.array([pdf_1,pdf_2,pdf_3,pdf_4,pdf_5])
# compute PDFs conflation
conflated_pdf = conflate(pdfs)
# ..............................
# >>> ========== <<< .
# plot !!!_NORMALIZED_!!! PDFs .
# >>> ========== <<< .
# ..............................
for i,pdf in enumerate(pdfs):
plt.plot(x,pdf/pdf.sum(),c=f'C{i}',label=f'PDF {i+1}',lw=1)
# normalize =============
# ^^^^^^^^^^^^^
# PDFs really do definitely need
# to be normalized,i.e. they have
# to sum to 1,because the cumulative
# probability needs to be 1 (100%)
# Plot conflated PDF
# NOTE: we don't need to normalize
# because it is already normalized
plt.plot(x,conflated_pdf,'k--',label='Conflated PDF',lw=2)
# Plot options here
plt.xlabel("x")
plt.ylabel("probability density")
plt.title("Conflated PDF")
plt.legend(loc='upper left')
plt.show()
