标签: poisson

我知道Knuth用于生成随机泊松分布数的算法(在Java下面)但是我如何将其转换为generateEvent()随机调用方法？

int poissonRandomNumber(int lambda) {
    double L = Math.exp(-lambda);
    int k = 0;
    double p = 1;
    do {
        k = k + 1;
        double u = Math.random();
        p = p * u;
    } while (p > L);
    return k - 1;
}

我需要模拟泊松等待时间.我发现了许多模拟到达次数的例子,但是我需要模拟一次到达的等待时间,给定平均等待时间.

我一直在寻找这样的代码:

public int getPoisson(double lambda) 
{   
    double L = Math.exp(-lambda);   
    double p = 1.0;   
    int k = 0;   

    do 
    {    
        k++;     
        p *= rand.nextDouble(); 
        p *= Math.random(); 
    } while (p > L);   

    return k - 1; 
} 

但这是到达人数,而不是到达时间.

效率优于精确度,更多是因为功耗比时间.我正在使用的语言是Java,如果算法只使用Random类中可用的方法,那将是最好的,但这不是必需的.

我正在尝试确定我可以准确地安排在C/C++中执行任务的粒度.目前我可以可靠地安排任务每5微秒发生一次,但我试图看看是否可以进一步降低这一点.

关于如何实现这一目标的任何建议/如果可能的话,将不胜感激.

因为我知道计时器粒度通常可以依赖于操作系统:我目前在Linux上运行,但如果时序粒度更好,则会使用Windows(虽然我不相信它,基于我在QueryPerformanceCounter中找到的)

我在裸机上执行所有测量(没有VM)./proc/timer_info确认我的CPU的纳秒定时器分辨率(但我知道这不会转换为纳秒警报分辨率)

当前

我现在的代码可以在这里找到

目前,我能够每5微秒(5000纳秒)执行一次请求,迟到的次数少于1%.当迟到确实发生时,它们通常仅落后一个周期(5000纳秒).

我现在正做三件事

将进程设置为实时优先级(这里有一些由@ Spudd86指出)

struct sched_param schedparm;
memset(&schedparm, 0, sizeof(schedparm));
schedparm.sched_priority = 99; // highest rt priority
sched_setscheduler(0, SCHED_FIFO, &schedparm);

最大限度地减少计时器松弛

prctl(PR_SET_TIMERSLACK, 1);

使用timerfds(2.6 Linux内核的一部分)

int timerfd = timerfd_create(CLOCK_MONOTONIC,0);
struct itimerspec timspec;
bzero(&timspec, sizeof(timspec));
timspec.it_interval.tv_sec = 0;
timspec.it_interval.tv_nsec = nanosecondInterval;
timspec.it_value.tv_sec = 0;
timspec.it_value.tv_nsec = 1;

timerfd_settime(timerfd, 0, &timspec, 0);

可能的改进

1. 将处理器专用于此过程？
2. 使用非阻塞timerfd,这样我就可以创建一个紧密循环,而不是阻塞(紧密循环会浪费更多CPU,但也可能更快响应警报)
3. 使用外部嵌入式设备进行触发(无法想象为什么会更好)

为什么

我目前正在为基准测试引擎创建工作负载生成器.工作负载生成器使用泊松过程模拟到达率(X请求/秒等).从泊松过程中,我可以确定必须从基准测试引擎发出请求的相对时间.

因此,例如,在每秒10次请求时,我们可能会发出以下请求:t = 0.02,0.04,0.05,0.056,0.09秒

这些请求需要提前安排然后执行.随着每秒请求数的增加,调度这些请求所需的粒度也会增加(每秒数千个请求需要亚毫秒精度).结果,我试图弄清楚如何进一步扩展这个系统.

我正在尝试使用XGBoost来模拟从不等长曝光时间段生成的数据的声明频率,但是无法使模型正确处理曝光.我通常会通过将log(曝光)设置为偏移量来实现此目的 - 您是否可以在XGBoost中执行此操作？

(这里发布了一个类似的问题:xgboost,偏移曝光？)

为了说明这个问题,下面的R代码使用以下字段生成一些数据:

• x1,x2 - 因子(0或1)
• 暴露 - 观察数据的政策期限
• 频率 - 每单位曝光的平均索赔数
• 声明 - 观察到的声明数量〜泊松(频率*暴露)

目标是使用x1和x2预测频率 - 真实模型是:如果x1 = x2 = 1则频率= 2,否则频率= 1.

曝光不能用于预测频率,因为在政策开始时不知道.我们可以使用它的唯一方法是:预期的索赔数量=频率*曝光率.

代码尝试使用XGBoost通过以下方式预测:

1. 将曝光设置为模型矩阵中的权重
2. 将日志(曝光)设置为偏移量

在这些下面,我已经展示了如何处理树(rpart)或gbm的情况.

set.seed(1)
size<-10000
d <- data.frame(
  x1 = sample(c(0,1),size,replace=T,prob=c(0.5,0.5)),
  x2 = sample(c(0,1),size,replace=T,prob=c(0.5,0.5)),
  exposure = runif(size, 1, 10)*0.3
)
d$frequency <- 2^(d$x1==1 & d$x2==1)
d$claims <- rpois(size, lambda = d$frequency * d$exposure)

#### Try to fit using XGBoost
require(xgboost)
param0 <- list(
  "objective"  = "count:poisson"
  , "eval_metric" = "logloss" …

我正在为 CPUE 数据运行零膨胀模型。该数据有零通货膨胀的证据，我已通过 Vuong 测试（在下面的代码中）确认了这一点。根据 AIC 的说法，完整模型 (zint) 优于零模型。我现在想要：

1. 检查完整模型的残差以确定模型拟合（由于缺乏来自同事、互联网和 R 书籍的信息而遇到麻烦）
2. 如果模型拟合看起来不错，则可视化模型的输出（使用偏移变量时如何制定实际参数值的方程）

我向该部门的几位统计学家寻求帮助（他们以前从未这样做过，并将我发送到相同的谷歌搜索网站），向统计部门本身（每个人都太忙）以及 stackoverflow feed 寻求帮助。

我很欣赏书籍的代码或指南（在线免费提供），其中包含使用偏移变量时处理可视化 ZIP 和模型拟合的代码。

 yc=read.csv("CPUE_ycs_trawl_withcobb_BLS.csv",header=TRUE)
 yc=yc[which(yc$countinyear<150),]
 yc$fyear=as.factor(yc$year_cap)
 yc$flocation=as.factor(yc$location)
 hist(yc$countinyear,20)
 yc$logoffset=log(yc$numtrawlyr)

 ###Run Zero-inflated poisson with offset for CPUE####
 null <- formula(yc$countinyear ~ 1| 1)
 znull <- zeroinfl(null, offset=logoffset,dist = "poisson",link = "logit",
 data = yc)

 int <- formula(yc$countinyear ~ assnage * spawncob| assnage * spawncob)
 zint <- zeroinfl(int, offset=logoffset,dist = "poisson",link = "logit", data  
 = yc)
 AIC(znull,zint)

  g1=glm(countinyear ~ assnage * spawncob,
  offset=logoffset,data=yc,family=poisson)
  summary(g1)

 ####Vuong …

我想构建一个随机森林回归器来模拟计数数据（泊松分布）。默认的“mse”损失函数不适合这个问题。有没有办法定义自定义损失函数并将其传递给 Python 中的随机森林回归器（Sklearn 等）？

是否有任何实现可以在任何包中拟合 Python 中的计数数据？

我正在使用 scipy stats 生成一个随机数。我使用了泊松分布。下面是一个例子：

import scipy.stats as sct

A =2.5
Pos = sct.poisson.rvs(A,size = 20)

当我打印 Pos 时，我得到以下数字：

array([1, 3, 2, 3, 1, 2, 1, 2, 2, 3, 6, 0, 0, 4, 0, 1, 1, 3, 1, 5])

从数组中可以看到生成了一些数字，例如6。

我想要限制最大数字（假设是 5），即使用 sct.poisson.rvs 生成的任何随机数都应该等于或小于 5，

我如何调整我的代码来实现它。顺便说一句，我在 Pandas Dataframe 中使用它。

我的问题是以最有效的方式提取N泊松随机值(RV),每个值具有不同的平均值/速率Lam.基本上是size(RV) == size(Lam).

这是一个天真(非常慢)的实现:

import numpy as NP

def multi_rate_poisson(Lam):
    rv = NP.zeros(NP.size(Lam))
    for i,lam in enumerate(Lam):
        rv[i] = NP.random.poisson(lam=lam, size=1)
    return rv

在我的笔记本电脑上,1e6样品给出:

Lam = NP.random.rand(1e6) + 1
timeit multi_poisson(Lam)
1 loops, best of 3: 4.82 s per loop

有可能从中改进吗？

我有一个问题,我想使用概率分布生成1到5之间的一组随机整数值.

Poisson和Inverse Gamma是两个分布,它们显示了我所发现的特征(多数均值,更低的数字).

我正在寻找使用Apache Commons Math,但我不知道如何使用可用的发行版生成我想要的数字.

我正在研究“使用指数平滑进行预测”。我被困在练习 16.4 的部分：

该数据集partx包含汽车零件的月销售历史记录。应用局部泊松模型。应该通过最大化似然或最小化平方误差总和来估计参数。

局部泊松模型定义为：

在哪里 和

我有以下代码，但似乎卡住了。优化总是返回接近起始值的东西。

我是否正确拟合了局部泊松模型？

library(expsmooth)
data("partx")
S <- function(x) {
  a <- x[1]
  if(a < 0 | a > 1)
    return(Inf)
  n <- length(partx)
  lambda <- numeric(n+1)
  error <- numeric(n)
  lambda[1] <- x[2]
  for(i in 1:n) {
    error[i] <- partx[i]-rpois(1,lambda[i])
    lambda[i+1] <-   (1-a)*lambda[i] + a*partx[i]
  }
  return(sum(error^2))
}

# returns a = 0.5153971 and lambda = 5.9282414
op1 <- optim(c(0.5,5L),S, control = list(trace = 1))
# returns a = 0.5999655 and lambda = …