总加权瑞利分布

编程问答 2022-05-01

问题描述

我根据帕累托分布函数用 nrow=1000（个体）和 ncol=100（天）模拟数据的步长：

set.seed(10)
sim_data <- replicate(100,VGAM::rpareto(1000,shape=7,scale=500))
sim_data <- as.data.frame(sim_data)
set.seed(10)
sim_data[,1:50] <- sim_data[50]*(-1) ##assign directionality
sim_data_directions <- as.data.frame(sim_data)
##randomize columns
sim_data <- sim_data_directions[,sample(ncol(sim_data_directions))] 
x <- c(1:31) ## for variable 'days' to associate to each step length
set.seed(10)
df_sim <- cbind(sim_data,t(apply(sim_data,1,function(x) {
  i1 <- sample(seq_along(x),1)
  out <- sum(sample(x,i1))
  c(days = i1,step_lengths = out)}
))) ## create step lengths
##adding weights for each level in variable days
df_sim <- as.data.frame(dplyr::add_count(df_sim,days))

使用此数据集 df_sim，模拟步长值、与每个步长相关的时间（以天为单位）和权重（每个时间变量的值数（以天为单位），我想总结分布，使用瑞利分布函数，其中每个天数水平的分布都被加权，如下所示：

rayleigh_distr <- sum(n*function (x) x*exp(-1*(x/2*sigma)^2)/sigma^2)

其中 n 是权重。 我如何根据权重总结每天的分布？

解决方法

首先，您的 Rayleigh PDF 中似乎有错误。应该是：

x*exp(-(x/sigma)^2/2)/sigma^2

看起来您希望 rayleigh_distr 作为函数返回具有权重 n 和尺度参数 sigma 的瑞利 mixture distribution 的 PDF。如果是这样，那就是（对于标量 x）：

rayleigh_distr <- function(x,sigma,n) sum(n*x*exp(-(x/sigma)^2/2)/sigma^2)

其中 sigma 和 n 是等长的向量。

如果您想为 x 传递一个向量（不一定与 sigma 和 n 的长度相同），那么这应该可行：

library(Rfast)
rayleigh_distr <- function(x,n) colsums(n*eachrow(exp(-(outer(1/sigma,x))^2/2),x)/sigma^2)

它将返回一个与 x 长度相同的向量。

更新：

CDF 将是：

rayleigh_cdf <- function(x,n) colSums(n*(1 - exp(-(outer(1/sigma,x))^2/2)))
distributiondistributionrrweighted