Как связать столбцы результатов моделирования Монте-Карло для логнормального распределения? - PullRequest
Новая
       106

Как связать столбцы результатов моделирования Монте-Карло для логнормального распределения?

0 голосов
/ 11 июля 2020

Я хочу выполнить логарифмическое нормальное моделирование Монте-Карло для 12 итераций с изменением расположения и формы распределения для каждой итерации. В качестве выходных данных следующего кода я хочу, чтобы каждая симуляция имела разные столбцы в виде фрейма данных. Например, если я хочу выполнить 10 000 симуляций для каждого набора местоположения и формы, результирующий фрейм данных будет иметь 12 столбцов и 10000 строк. Я написал код в двух частях, где сначала определил функцию для lapply, а затем попытался выполнить моделирование, при этом результаты каждой итерации сохраняются в другом столбце. Код части 1: 

myfxn <- function(i,s,location,shape){
cvsq=1
  i <-seq(5,60,5)
  s <-i*cvsq
  location <- log(i^2 / sqrt(s^2 + i^2))
  shape <- sqrt(log(1 + (s^2 / i^2)))
}

Код части 2: 

DF1 <- do.call(cbind, lapply(seq(5,60,5), function(myfxn) setNames(data.frame(rlnorm(n=10000, location, shape)), i)))
head(DF1)

Проблема в том, что я продолжаю получать ошибки с определенной функцией.

1 Ответ

0 голосов
/ 12 июля 2020

Кажется, вы запутались в объявлении вашей функции. Ваша функция выглядит так, как будто она должна принимать 4 параметра (i, s, местоположение и форма). Однако внутри функции эти входы никогда не используются - все они создаются с нуля в соответствии со значением cvsq. Поэтому они являются избыточными в качестве входных данных для функции.

Внутри вашего lapply вы передаете значения seq(5,60,5) вашей функции, но неясно, для чего они нужны, если переменные i, s , расположение и форма вычисляются заново каждый раз при вызове функции. Как вы хотите, чтобы ваша функция делала с этими числами? Если это так, то единственная переменная, которую имеет смысл передавать в качестве одного параметра, - это cvsq, поскольку изменение cvsq изменяет значения всех других переменных. Если это неверно, то просто непонятно, что вы пытаетесь сделать.

Если я на правильном пути, ваш код будет выглядеть так: 

myfxn <- function(cvsq)
{
  i <-seq(5,60,5)
  s <-i*cvsq
  location <- log(i^2 / sqrt(s^2 + i^2))
  shape <- sqrt(log(1 + (s^2 / i^2)))
  rlnorm(n = 10000, location, shape)
}

my_cvsq <- seq(5, 60, 5)
df <- as.data.frame(do.call(cbind, lapply(my_cvsq, myfxn)))

head(df)
#>           V1        V2         V3        V4         V5        V6           V7
#> 1  0.1786204 1.2587228 0.04588056 1.1925041 0.01985825 0.6924621 0.0008467456
#> 2 28.6147331 0.1243106 0.08306054 0.5057983 4.65819135 3.0851280 0.0429381886
#> 3  0.6666314 0.5546972 2.83435439 0.6512173 0.04994503 6.4365528 0.3779170170
#> 4  7.7163909 4.5954509 0.26637679 0.4085584 3.31589074 0.4894957 1.5195320720
#> 5  1.0931151 0.0575365 0.80689457 0.5033203 1.35200408 0.1427171 2.4058576253
#> 6  0.2573164 7.6190362 1.01222386 1.8598690 5.91525940 0.5600217 0.7858050404
#>            V8         V9        V10         V11       V12
#> 1 1.263350450 0.09485695 0.02615021 0.331885430 2.9109114
#> 2 0.004408080 0.14120821 0.05990997 0.053248650 0.1013357
#> 3 0.888602537 1.32448399 1.17625454 0.766153864 0.0196438
#> 4 5.254249758 0.01337194 0.08463637 2.996341976 0.8942556
#> 5 0.004375628 0.05052838 0.51024824 0.533356862 1.6166927
#> 6 0.985499853 0.09939285 0.86272478 0.001223073 1.0652496

И мы можем увидеть, как выглядят ваши 12 образцов, построив график их плотности: 

d <- density(df[,1], from = 0, to = 10)
plot(d$x, d$y, type = "l", xlim = c(0, 10), ylim = c(0, 0.7))
for(i in 2:12){
  d <- density(df[,i], from = 0, to = 10)
  lines(d$x, d$y, col = i)
}

Created on 2020-07-11 by the пакет реплекс (v0.3.0)

...