python之如何在 python 中权衡具有高斯分布的 2 个变量的函数

最近几天我一直在处理这个问题，但我还看不出问题出在哪里。

我正在尝试用具有特定平均值(R0) 和偏差 (sigma)。这是必需的，因为在实验分析时，理论函数 f(q) 在其 r 变量中具有一定的分散性。因此，我们使用概率密度函数来衡量我们在 r 变量中的功能。

我包含了有效的代码，但没有给出预期的结果(加权曲线应该随着多分散性的增长而变得更平滑(更高的 sigma)，如下所示。如您所愿看，我将 f(r,q)*g(r) 这两个函数的卷积从 r = 0 集成到 r = +inf.

绘制结果以将称重结果与简单函数进行比较:

from scipy.integrate import quad, quadrature 
import numpy as np 
import math as m 
import matplotlib.pyplot as plt  
 
#function weighted with a probability density function (gaussian) 
def integrand(r,q): 
    #gaussian function normalized 
    def gauss_nor(r): 
        #gaussian function 
        def gauss(r): 
            return m.exp(-((r-R0)**2)/(2*sigma**2)) 
        return (m.exp(-((r-R0)**2)/(2*sigma**2)))/(quad(gauss,0,np.inf)[0]) 
    #function f(r,q) 
    def f(r,q): 
        return 3*(np.sin(q*r)-q*r*np.cos(q*r))/((r*q)**3) 
    return gauss_nor(r)*f(r,q) 
 
#quadratic integration of the integrand (from 0 to +inf) 
#integrand is function*density_function (gauss)  
def function(q): 
    return quad(integrand, 0, np.inf, args=(q))[0] 
 
#parameters used in the function 
R0=20 
sigma=5 
 
#range to plot q 
q=np.arange(0.001,2.0,0.005) 
 
#vector where the result of the integral will be saved 
function_vec = np.vectorize(function) 
 
#vector for the squared power of the integral 
I=[] 
I=(function_vec(q))**2 
 
#function without density function 
I0=[] 
I0=(3*(np.sin(q*R0)-q*R0*np.cos(q*R0))/((R0*q)**3))**2 
 
#plot of weighted and non-weighted functions 
p1,=plt.plot(q,I,'b') 
p3,=plt.plot(q,I0,'r') 
plt.legend([p1,p3],('Weighted','No weighted')) 
plt.yscale('log') 
plt.xscale('log') 
plt.show() 

非常感谢。我已经遇到这个问题好几天了，但我还没有发现错误。

也许有人知道如何以更简单的方式用 PDF 权衡函数。

请您参考如下方法:

我简化了你的代码，输出结果和你的一样。我认为它已经很平滑了，log-log 图中有一些非常尖锐的峰值，只是因为曲线有零点。所以它在对数-对数图中并不平滑，但在正常的 X-Y 图中却很平滑。

import numpy as np 
 
def gauss(r): 
    return np.exp(-((r-R0)**2)/(2*sigma**2)) 
 
def f(r,q): 
    return 3*(np.sin(q*r)-q*r*np.cos(q*r))/((r*q)**3) 
 
R0=20 
sigma=5 
 
qm, rm = np.ogrid[0.001:2.0:0.005, 0.001:40:1000j] 
gr = gauss(rm) 
gr /= np.sum(gr) 
fm = f(rm, qm) 
fm *= gr 
 
plot(qm.ravel(), fm.sum(axis=1)**2) 
plt.yscale('log') 
plt.xscale('log')