образец логнормального распределения с точным средним и SD

Question

Я подготовил вектор путем выборки логарифмически нормального распределения путем (методом проб и ошибок) установки параметров для mean и sd, поэтому rlnorm() возвращает ровно среднее значение 20 и sd 6 (до 3 десятичных знаков) для любого указанного случайного значения set.seed() согласно следующему примеру ...

# 10,000 samples from log-normal distribution

set.seed(7)
HcT <- rlnorm(n = 10000, log(19.147), log(1.33832))

# Report mean and sd

paste('The mean of HcT is',round(mean(HcT),3),'and the SD is',round(sd(HcT),3)) 

[1] "The mean of HcT is 20 and the SD is 6"

Однако вместо проб и ошибок я хотел бы «стремиться к цели» по двум параметрам. Существует несколько примеров переполнения стека для поиска цели по одному значению , но я не уверен, какую функцию или пакет применить для случая с двумя параметрами (среднее значение и SD в приведенном выше примере).

Answer 1

Должно работать ОК, чтобы минимизировать сумму квадратов отклонений от целевых значений.В этом подходе есть подводные камни (см., Например, Числовые рецепты от Press et al.), Но он должен быть в порядке для простых задач.Следующий код появляется для получения правильных ответов для вашего случая:

f <- function(p,seed=7,target=c(20,6)) {
    mu <- log(p[1])
    sd <- log(p[2])
    set.seed(seed)
    r <- rlnorm(1e4,mu,sd)
    sum((c(mean(r),sd(r))-target)^2)
}

Выбор некоторых нелепых начальных значений ({15,2}):

optim(par=c(15,2), fn=f)

---

На основеПри ответе @ Коула я бы подумал, что это сработает отлично: нарисуйте нормальные отклонения, преобразуйте их так, чтобы они имели среднее значение и sd точно , равные значениям логарифмической шкалы, а затем возвели в степень.Но это работает только в среднем или асимптотически (т. Е. Большая выборка сходится к желаемому среднему значению), не совсем для конечных выборок.Точно не продумал, почему это так.

rlnorm_exact <- function(n, m, sd) {
    m2 <-  log(m^2 / sqrt(sd^2 + m^2))
    sd2 <- sqrt(log(1 + (sd^2 / m^2)))
    r <- c(scale(rnorm(n)))
    return(exp(sd2*r+m2))
}

Answer 2

Я тоже сталкивался с этой проблемой и раньше, и ссылка ниже прояснила меня.rlnorm() не просто использует логарифм среднего арифметического и стандартного отклонения.Вместо этого функция ожидает mu и sigma, специфичные для логнормального распределения.

К счастью, люди по этой ссылке вывели формулы для нас, чтобы мы преобразовали их в логнормальные распределения.

Я собираюсь сделать это менее привлекательным, чтобы люди обращались кссылка выше, как они решили это:

m <- 20
s <- 6

data_set <- rlnorm(n=1000000, 
                   meanlog=log(m^2 / sqrt(s^2 + m^2)), 
                   sdlog=sqrt(log(1 + (s^2 / m^2))))

mean(data_set)
sd(data_set)

Редактировать: измененная переменная с sd на s, потому что sd() также является функцией ...