- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
我正在尝试重现 http://jheusser.github.io/2013/09/08/hawkes.html 中的作品在 python 中,除了不同的数据。我已经编写了代码来模拟泊松过程以及他们描述的霍克斯过程。
为了执行 Hawkes 模型 MLE,我将对数似然函数定义为
def loglikelihood(params, data):
(mu, alpha, beta) = params
tlist = np.array(data)
r = np.zeros(len(tlist))
for i in xrange(1,len(tlist)):
r[i] = math.exp(-beta*(tlist[i]-tlist[i-1]))*(1+r[i-1])
loglik = -tlist[-1]*mu
loglik = loglik+alpha/beta*sum(np.exp(-beta*(tlist[-1]-tlist))-1)
loglik = loglik+np.sum(np.log(mu+alpha*r))
return -loglik
使用一些虚拟数据,我们可以计算 Hawkes 过程的 MLE
atimes=[58.98353497, 59.28420225, 59.71571013, 60.06750179, 61.24794134,
61.70692463, 61.73611983, 62.28593814, 62.51691723, 63.17370423
,63.20125152, 65.34092403, 214.24934446, 217.0390236, 312.18830525,
319.38385604, 320.31758188, 323.50201334, 323.76801537, 323.9417007]
res = minimize(loglikelihood, (0.01, 0.1,0.1),method='Nelder-Mead',args = (atimes,))
print res
但是,我不知道如何在 python 中执行以下操作。
最佳答案
好的,您可能希望做的第一件事就是绘制数据。为了简单起见,我转载了 this figure因为它只发生了 8 个事件,所以很容易看到系统的行为。以下代码:
import numpy as np
import math, matplotlib
import matplotlib.pyplot
import matplotlib.lines
mu = 0.1 # Parameter values as found in the article http://jheusser.github.io/2013/09/08/hawkes.html Hawkes Process section.
alpha = 1.0
beta = 0.5
EventTimes = np.array([0.7, 1.2, 2.0, 3.8, 7.1, 8.2, 8.9, 9.0])
" Compute conditional intensities for all times using the Hawkes process. "
timesOfInterest = np.linspace(0.0, 10.0, 100) # Times where the intensity will be sampled.
conditionalIntensities = [] # Conditional intensity for every epoch of interest.
for t in timesOfInterest:
conditionalIntensities.append( mu + np.array( [alpha*math.exp(-beta*(t-ti)) if t > ti else 0.0 for ti in EventTimes] ).sum() ) # Find the contributions of all preceding events to the overall chance of another one occurring. All events that occur after t have no contribution.
" Plot the conditional intensity time history. "
fig = matplotlib.pyplot.figure()
ax = fig.gca()
labelsFontSize = 16
ticksFontSize = 14
fig.suptitle(r"$Conditional\ intensity\ VS\ time$", fontsize=20)
ax.grid(True)
ax.set_xlabel(r'$Time$',fontsize=labelsFontSize)
ax.set_ylabel(r'$\lambda$',fontsize=labelsFontSize)
matplotlib.rc('xtick', labelsize=ticksFontSize)
matplotlib.rc('ytick', labelsize=ticksFontSize)
eventsScatter = ax.scatter(EventTimes,np.ones(len(EventTimes))) # Just to indicate where the events took place.
ax.plot(timesOfInterest, conditionalIntensities, color='red', linestyle='solid', marker=None, markerfacecolor='blue', markersize=12)
fittedPlot = matplotlib.lines.Line2D([],[],color='red', linestyle='solid', marker=None, markerfacecolor='blue', markersize=12)
fig.legend([fittedPlot, eventsScatter], [r'$Conditional\ intensity\ computed\ from\ events$', r'$Events$'])
matplotlib.pyplot.show()
相当准确地重现了这个数字,尽管我选择的事件时期有些武断:
这也可以应用于 example set of data 的集合5000 笔交易,方法是对数据进行分箱并将每个分箱视为一个事件。然而,现在发生的情况是,每个事件的权重都略有不同,因为每个容器中发生的交易数量不同。the article中也提到了这一点在将比特币交易到达与霍克斯过程相匹配部分提出了克服这个问题的方法:与原始数据集的唯一区别是我为所有共享一个交易的交易添加了一个随机毫秒时间戳与另一笔交易的时间戳。这是必需的,因为模型需要区分每笔交易(即每笔交易必须具有唯一的时间戳)。
这包含在以下代码中:
import numpy as np
import math, matplotlib, pandas
import scipy.optimize
import matplotlib.pyplot
import matplotlib.lines
" Read example trades' data. "
all_trades = pandas.read_csv('all_trades.csv', parse_dates=[0], index_col=0) # All trades' data.
all_counts = pandas.DataFrame({'counts': np.ones(len(all_trades))}, index=all_trades.index) # Only the count of the trades is really important.
empirical_1min = all_counts.resample('1min', how='sum') # Bin the data so find the number of trades in 1 minute intervals.
baseEventTimes = np.array( range(len(empirical_1min.values)), dtype=np.float64) # Dummy times when the events take place, don't care too much about actual epochs where the bins are placed - this could be scaled to days since epoch, second since epoch and any other measure of time.
eventTimes = [] # With the event batches split into separate events.
for i in range(len(empirical_1min.values)): # Deal with many events occurring at the same time - need to distinguish between them by splitting each batch of events into distinct events taking place at almost the same time.
if not np.isnan(empirical_1min.values[i]):
for j in range(empirical_1min.values[i]):
eventTimes.append(baseEventTimes[i]+0.000001*(j+1)) # For every event that occurrs at this epoch enter a dummy event very close to it in time that will increase the conditional intensity.
eventTimes = np.array( eventTimes, dtype=np.float64 ) # Change to array for ease of operations.
" Find a fit for alpha, beta, and mu that minimises loglikelihood for the input data. "
#res = scipy.optimize.minimize(loglikelihood, (0.01, 0.1,0.1), method='Nelder-Mead', args = (eventTimes,))
#(mu, alpha, beta) = res.x
mu = 0.07 # Parameter values as found in the article.
alpha = 1.18
beta = 1.79
" Compute conditional intensities for all epochs using the Hawkes process - add more points to see how the effect of individual events decays over time. "
conditionalIntensitiesPlotting = [] # Conditional intensity for every epoch of interest.
timesOfInterest = np.linspace(eventTimes.min(), eventTimes.max(), eventTimes.size*10) # Times where the intensity will be sampled. Sample at much higher frequency than the events occur at.
for t in timesOfInterest:
conditionalIntensitiesPlotting.append( mu + np.array( [alpha*math.exp(-beta*(t-ti)) if t > ti else 0.0 for ti in eventTimes] ).sum() ) # Find the contributions of all preceding events to the overall chance of another one occurring. All events that occur after time of interest t have no contribution.
" Compute conditional intensities at the same epochs as the empirical data are known. "
conditionalIntensities=[] # This will be used in the QQ plot later, has to have the same size as the empirical data.
for t in np.linspace(eventTimes.min(), eventTimes.max(), eventTimes.size):
conditionalIntensities.append( mu + np.array( [alpha*math.exp(-beta*(t-ti)) if t > ti else 0.0 for ti in eventTimes] ).sum() ) # Use eventTimes here as well to feel the influence of all the events that happen at the same time.
" Plot the empirical and fitted datasets. "
fig = matplotlib.pyplot.figure()
ax = fig.gca()
labelsFontSize = 16
ticksFontSize = 14
fig.suptitle(r"$Conditional\ intensity\ VS\ time$", fontsize=20)
ax.grid(True)
ax.set_xlabel(r'$Time$',fontsize=labelsFontSize)
ax.set_ylabel(r'$\lambda$',fontsize=labelsFontSize)
matplotlib.rc('xtick', labelsize=ticksFontSize)
matplotlib.rc('ytick', labelsize=ticksFontSize)
# Plot the empirical binned data.
ax.plot(baseEventTimes,empirical_1min.values, color='blue', linestyle='solid', marker=None, markerfacecolor='blue', markersize=12)
empiricalPlot = matplotlib.lines.Line2D([],[],color='blue', linestyle='solid', marker=None, markerfacecolor='blue', markersize=12)
# And the fit obtained using the Hawkes function.
ax.plot(timesOfInterest, conditionalIntensitiesPlotting, color='red', linestyle='solid', marker=None, markerfacecolor='blue', markersize=12)
fittedPlot = matplotlib.lines.Line2D([],[],color='red', linestyle='solid', marker=None, markerfacecolor='blue', markersize=12)
fig.legend([fittedPlot, empiricalPlot], [r'$Fitted\ data$', r'$Empirical\ data$'])
matplotlib.pyplot.show()
这会生成以下与绘图的拟合: 一切看起来都不错,但是,当您查看细节时,您会发现通过简单地获取交易数量的一个向量并减去拟合的向量来计算残差是行不通的,因为它们具有不同的长度: 然而,可以在与记录经验数据时相同的时期提取强度,然后计算残差。这使您能够找到经验数据和拟合数据的分位数并将它们相互绘制,从而生成 QQ 图:
""" GENERATE THE QQ PLOT. """
" Process the data and compute the quantiles. "
orderStatistics=[]; orderStatistics2=[];
for i in range( empirical_1min.values.size ): # Make sure all the NANs are filtered out and both arrays have the same size.
if not np.isnan( empirical_1min.values[i] ):
orderStatistics.append(empirical_1min.values[i])
orderStatistics2.append(conditionalIntensities[i])
orderStatistics = np.array(orderStatistics); orderStatistics2 = np.array(orderStatistics2);
orderStatistics.sort(axis=0) # Need to sort data in ascending order to make a QQ plot. orderStatistics is a column vector.
orderStatistics2.sort()
smapleQuantiles=np.zeros( orderStatistics.size ) # Quantiles of the empirical data.
smapleQuantiles2=np.zeros( orderStatistics2.size ) # Quantiles of the data fitted using the Hawkes process.
for i in range( orderStatistics.size ):
temp = int( 100*(i-0.5)/float(smapleQuantiles.size) ) # (i-0.5)/float(smapleQuantiles.size) th quantile. COnvert to % as expected by the numpy function.
if temp<0.0:
temp=0.0 # Avoid having -ve percentiles.
smapleQuantiles[i] = np.percentile(orderStatistics, temp)
smapleQuantiles2[i] = np.percentile(orderStatistics2, temp)
" Make the quantile plot of empirical data first. "
fig2 = matplotlib.pyplot.figure()
ax2 = fig2.gca(aspect="equal")
fig2.suptitle(r"$Quantile\ plot$", fontsize=20)
ax2.grid(True)
ax2.set_xlabel(r'$Sample\ fraction\ (\%)$',fontsize=labelsFontSize)
ax2.set_ylabel(r'$Observations$',fontsize=labelsFontSize)
matplotlib.rc('xtick', labelsize=ticksFontSize)
matplotlib.rc('ytick', labelsize=ticksFontSize)
distScatter = ax2.scatter(smapleQuantiles, orderStatistics, c='blue', marker='o') # If these are close to the straight line with slope line these points come from a normal distribution.
ax2.plot(smapleQuantiles, smapleQuantiles, color='red', linestyle='solid', marker=None, markerfacecolor='red', markersize=12)
normalDistPlot = matplotlib.lines.Line2D([],[],color='red', linestyle='solid', marker=None, markerfacecolor='red', markersize=12)
fig2.legend([normalDistPlot, distScatter], [r'$Normal\ distribution$', r'$Empirical\ data$'])
matplotlib.pyplot.show()
" Make a QQ plot. "
fig3 = matplotlib.pyplot.figure()
ax3 = fig3.gca(aspect="equal")
fig3.suptitle(r"$Quantile\ -\ Quantile\ plot$", fontsize=20)
ax3.grid(True)
ax3.set_xlabel(r'$Empirical\ data$',fontsize=labelsFontSize)
ax3.set_ylabel(r'$Data\ fitted\ with\ Hawkes\ distribution$',fontsize=labelsFontSize)
matplotlib.rc('xtick', labelsize=ticksFontSize)
matplotlib.rc('ytick', labelsize=ticksFontSize)
distributionScatter = ax3.scatter(smapleQuantiles, smapleQuantiles2, c='blue', marker='x') # If these are close to the straight line with slope line these points come from a normal distribution.
ax3.plot(smapleQuantiles, smapleQuantiles, color='red', linestyle='solid', marker=None, markerfacecolor='red', markersize=12)
normalDistPlot2 = matplotlib.lines.Line2D([],[],color='red', linestyle='solid', marker=None, markerfacecolor='red', markersize=12)
fig3.legend([normalDistPlot2, distributionScatter], [r'$Normal\ distribution$', r'$Comparison\ of\ datasets$'])
matplotlib.pyplot.show()
这会生成以下图:
经验数据的分位数图并不完全是 the same as in the article ,我不确定为什么,因为我不擅长统计。但是,从编程的角度来看,这就是您完成所有这些工作的方式。
关于python - 如何计算python中点过程的残差,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/24785518/
如果我声明了类似的类型 type test(NSIZE) integer, len :: NSIZE real :: dummy(NSIZE) contains procedure,
我知道这是一个不太可能的事情,但是由于“选项私有(private)模块”的限制,甚至更糟糕的“私有(private)子/函数”的限制,有谁知道是否有一种方法可以从 Excel 应用程序隐藏 VBA 过
我有两个表,property 和 component。 component.id_property = property.id。 我正在尝试创建一个过程,该过程对所选属性的组件进行计数,如果所选属性没
我有一份报告,它是在 SSRS 2005 中开发的,我正在使用存储过程从数据库中获取结果。报告输出的结果非常简单,如下图所示。 如果假设我正在寻找不同的成员 例如:- MemberID c108 c
我需要一个通用函数/过程,该函数/过程将根据提供的数据计算出我的淡入淡出时间和值,如下所示: 我将字节值保存在字节数组中:这些是起始值。然后,我在其他数组中存储了一些值:这些将是新值。然后我有时间要提
我想在界面的多个按钮上创建相同的操作。是否只能通过创建单独的操作监听器方法并调用执行操作的方法才可行,还是还有其他方法?是否可以将按钮放在一个组中并执行以下操作:- groupButton.setOn
我有以下情况: procedure Test; begin repeat TryAgain := FALSE; try // Code // Code if this an
我正在尝试执行以下操作;假设我在 Oracle 中创建了一个对象类型 create type test as object( name varchar2(12), member procedure p
问题: 如果可能的话,如何声明一个用于任何类型参数的函数 T其中 T 的唯一约束是它被定义为 1D array如 type T is array ( integer range <> ) of a_r
我正在尝试创建这个 mysql 过程来制作一个包含今年所有日期和所有时间的表(以一小时为间隔。) CREATE TABLE FECHAS ( created_at datetime ); CREA
所以, 我在这里面临一个问题,这让我发疯,我认为这是一个愚蠢的错误,所以我不是 MySQL 的新手,但它并不像我想象的那样工作。 尝试将此语句部署到 MySQL 后,我收到此错误: ERROR 106
我有一个架构,其中包含星球大战中的人物列表、他们出现的电影、他们访问的行星等。这是架构: CREATE DATABASE IF NOT EXISTS `starwarsFINAL` /*!40100
我一直在为一家慈善机构创建一款应用程序,允许家庭在节日期间注册接收礼物。数据库组织有多个表。下面列出了这些表(及其架构/创建语句): CREATE TABLE IF NOT EXISTS ValidD
正如上面标题所解释的,我正在尝试编写一个sql函数来按日期删除表而不删除系统表。我在此消息下方放置了一张图片,以便直观地解释我的问题。任何帮助将不胜感激!感谢您的时间! 最佳答案 您可以通过查询INF
DELIMITER $$ CREATE PROCEDURE INSERT_NONE_HISTORY_CHECKBOX() BEGIN DECLARE note_id bigint(20); F
是否可以编写一个存储过程或触发器,在特定时间在数据库内部自动执行,而无需来自应用程序的任何调用?如果是,那么任何人都可以给我一个例子或链接到一些我可以阅读如何做到这一点的资源。 最佳答案 查看 pgA
我需要创建一个过程:1)从表中的字段中选择一些文本并将其存储在变量中2) 更新相同的记录字段,仅添加 yyyymmdd 格式的日期以及过程中的附加文本输入...类似这样的... delimiter /
好的,这就是我想做的: 如果条目已存在(例如基于字段name),则只需返回其id 如果没有,请添加 这是我迄今为止所管理的(对于“如果不存在,则创建它”部分): INSERT INTO `object
以下是我编写的程序,用于找出每位客户每天购买的前 10 件商品。 这是我尝试过的第一个 PL/SQL 操作。它没有达到我预期的效果。 我使用的逻辑是接受开始日期、结束日期以及我对每个客户感兴趣的前“x
我正在尝试在MySQL中创建一个过程那insert week s(当年)发送至我的 week table 。但存在一个问题,因为在为下一行添加第一行后,我收到错误: number column can
我是一名优秀的程序员,十分优秀!