2

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

我正在尝试在具有特定平均值 ( ) 和偏差 ( )f(q,r)的高斯分布中使用 2 个变量对函数进行加权。这是必需的，因为理论函数在实验分析时其变量具有一定的分散性。因此，我们使用概率密度函数来衡量我们在变量中的函数。g(r)R0sigmaf(q)rr

我包含了有效的代码，但没有给出预期的结果（加权曲线应该随着多分散性的增长（更高sigma）更平滑，如下所示。如您所见，我整合了来自的 2 个函数的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 来衡量函数。

4

1 回答 1

1

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

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')

在此处输入图像描述

于 2013-04-02T11:27:49.750 回答