问题描述
我已经使用盒式积分生成了一些数据，基于它们是否位于直线之下。生成此数据的代码如下：
class Linear:
def __init__(self,xmin,xmax,intercept,slope):
self.xmin=xmin
self.xmax=xmax
self.intercept=intercept
self.slope=slope
self.mass=[]
def evaluate_linear(self,t):
return self.intercept + self.slope*t
def maxval(self):
return self.evaluate_linear(0)
def next_linear(self):
do_it=True
while(do_it==True):
x = np.random.uniform(self.xmin,self.xmax)
y1 = self.evaluate_linear(x)
maximum=self.evaluate_linear(0)
y2 = np.random.uniform(0,self.evaluate_linear(0))
if (y2 < y1):
filtered_x = x
self.mass.append(filtered_x)
return filtered_x
def integrate_linear(self):
func = lambda t: self.intercept + self.slope*t
return integrate.quad(func,self.xmin,self.xmax)
class Signal:
def __init__(self,mean,sigma):
self.xmin= xmin
self.xmax=xmax
self.signal_mean= mean
self.signal_sigma=sigma
self.mass=[]
def evaluate_gaussian(self,t):
return (1/np.sqrt(2*np.pi*self.signal_sigma**2))*np.exp(-((t-self.signal_mean)**2)/(2*self.signal_sigma**2))
def next_gaussian(self):
doLoop = True
while(doLoop==True):
x = np.random.uniform(self.xmin,self.xmax)
y1 = self.evaluate_gaussian(x)
y2 = np.random.uniform(0,1/np.sqrt(2*np.pi*self.signal_sigma**2))
if (y2<y1):
filtered_x = x
self.mass.append(filtered_x)
return filtered_x
def integrate_gaussian(self):
func = lambda t: (1/np.sqrt(2*np.pi*self.signal_sigma**2))*np.exp(-((t-self.signal_mean)**2)/(2*self.signal_sigma**2))
return integrate.quad(func,self.xmax)
class SignalwithBackground:
def __init__(self,sigma,sig_fraction,slope):
self.signal_mean= mean
self.signal_sigma=sigma
self.sig_fraction = sig_fraction
self.intercept=intercept
self.slope=slope
self.xmin= xmin
self.xmax= xmax
self.mass=[]
self.mass_sig=[]
self.mass_bgd=[]
def next(self):
q = np.random.uniform()
if (q < self.sig_fraction):
filtered_x = Signal(self.xmax,self.signal_mean,self.signal_sigma).next_gaussian()
self.mass_sig.append(filtered_x)
else:
filtered_x = Linear(self.xmin,self.xmax,self.intercept,self.slope).next_linear()
self.mass_bgd.append(filtered_x)
self.mass.append(filtered_x)
return filtered_x
def Plot(xmin,slope,NBINS,nevents_sig,nevents_bgd,plotting):
#xmin=0
#xmax=20
#intercept=20
#slope=-1
#mean=10
#sigma=0.5
#NBINS=100
sig_fraction = nevents_sig/(nevents_bgd + nevents_sig)
pdf = SignalwithBackground(xmax,slope)
for i in range( nevents_sig + nevents_bgd): pdf.next()
#print(data)
data= pdf.mass
sig_data = pdf.mass_sig
bgd_data = pdf.mass_bgd
if plotting==True:
myRange = (xmin,xmax)
fig,axs = plt.subplots (3,1,sharex='col')
axs[0].set_title("Signal distribution (" + str(len(sig_data )) + " entries)")
axs[1].set_title("Background distribution (" + str(len(bgd_data )) + " entries)")
axs[2].set_title("Total distribution (" + str(len(data)) + " entries)")
axs[2].set_xlabel('X')
axs[0].hist(sig_data,bins=NBINS,range=myRange)
axs[1].hist(bgd_data,bins = NBINS)
axs[2].hist(data,bins = NBINS)
fig.tight_layout ()
plt.savefig('Example1.pdf')
else:
return data,sig_data,bgd_data
然后，我希望拟合我的数据，以便使用scipy.optimize.minimise将对数似然的卡方等效值最大化，将我的期望值作为用于生成分布的线性函数的值，在每个垃圾箱的中点。当我运行代码时，没有出现这样的错误，但是即使更改了用于生成分布的线性函数，我希望最小化的值也没有改变，这使我相信自己做错了。的代码如下：
def log_likelihood(vars,nobserved):
data=[]
data2=np.linspace(20.2,0.2,100)
data2*vars[0] + vars[1]
for i in range(100):
if data2[i]>0:
data.append(-2*(data2[i]-nobserved[i]+nobserved[i]*np.log(nobserved[i]/data2[i])))
else:
data.append(0)
return sum(data)
data,signal,background = Plot(0,20,15,-1,10,0.5,100,150,10000,False)
vars=[-1,20]
result = scipy.optimize.minimize(log_likelihood,vars,args =(background))
print(result)
任何帮助将不胜感激！
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）