R에서 연속 분포 매개변수 추정을 위한 모의 어닐링

이것은 재미있는 것이었다!

일반 인구에서 설명하는 분포HbA1c를 찾으려고 노력했습니다. 몇 시간 동안 파헤친 후 내가 찾을 수 있는 최선은 NIH의 이 연구였습니다. 이 기사에서 우리에게 주어진 것은 몇 개의 분위수입니다.


5%
10%
25%
50%
75%
90%
95%


4.7
4.8
5.1
5.4
5.7
6.1
6.8


이 점에 연속 분포를 맞추려고 합니다. 모의 어닐링을 통해 문제 공간에서 전역 최적값을 검색할 수 있습니다. 연구 텍스트에서 데이터가 왜곡되어 있다는 것을 알고 있습니다(평균 = 5.6, 모드 = 5.3). 왜도를 허용하는 분포가 필요합니다. gamma distribution으로 갑시다.

가장 먼저 해야 할 일은 알려진 데이터 세트의 데이터 프레임을 만드는 것입니다.

df <- data.frame(quantile = c(0.05, 0.1, 0.25, 0.50, 0.75, 0.9, 0.95),
                 value = c(4.7, 4.8, 5.1, 5.4, 5.7, 6.1, 6.8))

다음으로 SA 알고리즘이 최소화할 함수를 작성해야 합니다.

fun <- function (x) {
    shape <- x[1]
    scale <- x[2]
    df$est <- qgamma(df$quantile, shape, scale)
    df$sq_err <- (df$est - df$value) ^ 2
    error <- sum(df$sq_err)
    error
}

모의 어닐링을 위해 GenSA 패키지를 사용할 것입니다.

library(GenSA)

lower_bound <- 1e-6
upper_bound <- 100

fit <- GenSA(par = c(1, 1), 
             fn = fun, 
             lower = c(lower_bound, lower_bound), 
             upper = c(upper_bound, upper_bound))

이는 모양 = 89.4837 및 규모 = 16.2514의 감마 분포 매개변수를 제공합니다. 이 피팅의 SSE는 0.1747입니다.

원본 데이터에 대해 플로팅했을 때 이것이 어떻게 보이는지 봅시다.

library(ggplot2)

df2 <- data.frame(x = seq(from = 3.0, to = 10.0, by = 0.01))
df2$y = pgamma(df2$x, fit$par[1], fit$par[2])

ggplot() +
    geom_line(aes(x, y), df2, color = "red") +
    geom_point(aes(value, quantile), df, size = 2, color = "blue")

자! 몇 개의 데이터 포인트와 분위수만 있는 경우 해당 분포를 추측할 수 있습니다. 여러 분포로 이것을 시도하고 결과 SSE를 비교할 수 있습니다.