如何修改从pdf生成随机数的功能更快运行
How to modify the function of generating random numbers from a pdf to run faster
下面是从pdf生成随机数的代码:
我修改了函数 rcmp(包 COMPoissonReg)的代码。
dcomp <- function(y,mu,v,z=NULL, max=100)
{
if (is.null(z)){
z=sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
}
log.ff <- v*y*log(mu)-v*lgamma(y) - log(z)
return(exp(log.ff))
}
rcomp <- function(n, mu, v, max=100)
{
if (length(mu) == 1) {
mu <- rep(mu, n)
}
if (length(v) == 1) {
v <- rep(v, n)
}
u <- runif(n)
y <- numeric(n)
z=sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
for (i in 1:n) {
px <- dcomp(y[i], mu[i], v[i],z=z[i], max = max)
while (px < u[i]) {
y[i] <- y[i] + 1
px <- px + dcomp(y[i], mu[i], v[i], z=z[i],max = max)
}
}
return(y)
}
但是,该函数花了很长时间来模拟随机变量,有没有办法将此代码修改为更快运行?
您的代码中存在一些错误,导致您的实施需要很长时间。
密度函数dcomp
修改如下
dcomp <- function(y,mu,v,z=NULL, max=100) {
if (is.null(z)){
z=sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
}
log.ff <- v*y*log(mu)-v*lgamma(y+1) - log(z)
return(exp(log.ff))
}
请注意,您需要将 lgamma
加 1 作为 gamma(x+1) = factorial(x)。
在生成随机变量的函数 rcomp
中,您在
中遇到了问题
sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
此行中的 sum
折叠矢量化。您需要更新它以获得具有 z
个性化值的适当向量。我在下面的代码中删除了这个预先计算,只在 dcomp
中进行计算,但是进行预先计算肯定会节省时间。
更新后的 rcomp
函数如下所示
rcomp <- function(n, mu, v, max=100)
{
if (length(mu) == 1) {
mu <- rep(mu, n)
}
if (length(v) == 1) {
v <- rep(v, n)
}
u <- runif(n)
y <- rep(0, n) # Have changed this line to force zeros as starting points
# z=sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
for (i in 1:n) {
px <- dcomp(y[i], mu[i], v[i], max = max) # Not using z
while (px < u[i]) {
y[i] <- y[i] + 1
px <- px + dcomp(y[i], mu[i], v[i],max = max) # Also not using z
}
}
return(y)
}
希望对您有所帮助!
下面是从pdf生成随机数的代码:
我修改了函数 rcmp(包 COMPoissonReg)的代码。
dcomp <- function(y,mu,v,z=NULL, max=100)
{
if (is.null(z)){
z=sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
}
log.ff <- v*y*log(mu)-v*lgamma(y) - log(z)
return(exp(log.ff))
}
rcomp <- function(n, mu, v, max=100)
{
if (length(mu) == 1) {
mu <- rep(mu, n)
}
if (length(v) == 1) {
v <- rep(v, n)
}
u <- runif(n)
y <- numeric(n)
z=sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
for (i in 1:n) {
px <- dcomp(y[i], mu[i], v[i],z=z[i], max = max)
while (px < u[i]) {
y[i] <- y[i] + 1
px <- px + dcomp(y[i], mu[i], v[i], z=z[i],max = max)
}
}
return(y)
}
但是,该函数花了很长时间来模拟随机变量,有没有办法将此代码修改为更快运行?
您的代码中存在一些错误,导致您的实施需要很长时间。
密度函数dcomp
修改如下
dcomp <- function(y,mu,v,z=NULL, max=100) {
if (is.null(z)){
z=sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
}
log.ff <- v*y*log(mu)-v*lgamma(y+1) - log(z)
return(exp(log.ff))
}
请注意,您需要将 lgamma
加 1 作为 gamma(x+1) = factorial(x)。
在生成随机变量的函数 rcomp
中,您在
sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
此行中的
sum
折叠矢量化。您需要更新它以获得具有z
个性化值的适当向量。我在下面的代码中删除了这个预先计算,只在dcomp
中进行计算,但是进行预先计算肯定会节省时间。
更新后的 rcomp
函数如下所示
rcomp <- function(n, mu, v, max=100)
{
if (length(mu) == 1) {
mu <- rep(mu, n)
}
if (length(v) == 1) {
v <- rep(v, n)
}
u <- runif(n)
y <- rep(0, n) # Have changed this line to force zeros as starting points
# z=sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))
for (i in 1:n) {
px <- dcomp(y[i], mu[i], v[i], max = max) # Not using z
while (px < u[i]) {
y[i] <- y[i] + 1
px <- px + dcomp(y[i], mu[i], v[i],max = max) # Also not using z
}
}
return(y)
}
希望对您有所帮助!