Расчет куртоза вручную против числовой интеграции - PullRequest
Новая
       30

Расчет куртоза вручную против числовой интеграции

1 голос
/ 04 ноября 2019

Я рассчитываю эксцесс функции плотности вероятности двумя разными методами.

У меня есть 2D-карты нормированного удельного углового момента j/jmean, дисперсии скорости s и массы m по радиусам r. Дисперсия скорости в данном пикселе i приводит к распределению j в этом пикселе, PDF_i, со средним j/jmean и стандартным отклонением r*s/jmean.

. взвешенное по массе среднее значение PDF_i во всех пикселях, PDF = sum(m*PDF_i)/sum(m).

Метод 1 включает в себя отбор выборки для каждого пикселя, взвешивание по массе, а затем вычисление эксцесса всего набора с помощью cumstats илипакет моментов.

Метод 2 включает вычисление эксцесса с помощью численного интегрирования: mu_4 = integrate((x-mu)^4 * PDF(x)dx).

Мой минимальный рабочий пример следующий: 

                                            # test methods for measuring kurtosis of a random distribution.
library(cumstats)
library(moments)
                                        # make a distribution of j, r, sigma, m
set.seed(1)
num <- 100
r <- runif(num)
rd <- 1.5
vflat <- 200
rflat <- 2
v <- vflat*(1-exp(-r/rflat))
s <- exp(-r/rd)
m <- exp(-r/rd)
j <- r*v
jmean <- sum(j*m)/sum(m)
j <- j/jmean
rs <- r*s/jmean

                                        # method 1:
                                        # in every spaxel, generate a PDF(j/jmean) with mean j and width r*s
res <- 300 #samples per spaxel
rho_j  <- rep(NA, res*length(j))
rho_m  <- rep(NA, res*length(m))
for (i in seq(length(j))) {
    rho_j[((i-1)*res+1):(i*res)] <-
        rnorm(res, mean=j[i], sd=(rs[i]))
    rho_m[((i-1)*res+1):(i*res)] <- rep(m[i], res)
}
                                        # apply a mass weighting
rho_jw <- rho_j*rho_m/(sum(rho_m)/num/res) #rho_j*rho_m gives the same results
                                        # use cumstats or moments to measure kurtosis of the whole thing.
b2cs <- cumstats::kurtosis(rho_jw)
b2mm <- moments::kurtosis(rho_jw)

print("res")
print(b2cs)
print(b2mm) # these are the same.

                                        # method 2:
                                        # apply Eq. 8 from paper
aint <- max(j+3*r*s) #integration limits +-A

integ <- function(x, mu, n) {
    (x - mu)^n*sum(m*dnorm(x, mean=j, sd=rs))/sum(m)
}
integrand <- function(x, mu, n) {sapply(x, integ, mu=mu, n=n)}

mu1 <- integrate(integrand, -aint, aint, mu=0, n=1)$value
sd <- sqrt(integrate(integrand, -aint, aint, mu=mu1, n=2)$value)
mu4 <- integrate(integrand, -aint, aint, mu=mu1, n=4)$value
eq8 <- mu4/sd^4

print("eq8")
print(eq8)

Я ожидаю, что два методадать тот же результат, но я получаю разные результаты (Res = 1,68; EQ8 = 2,21). Может показаться, что это вызвано моей реализацией, но я не вижу, в чем дело.

...