I get two similar graphs.
https://web.ncf.ca/nashjc/jfiles/Rplot-Labone-4095.pdf
https://web.ncf.ca/nashjc/jfiles/RplotLabone10K.pdf
Context:
R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 20.2
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
locale:
[1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C LC_TIME=en_CA.UTF-8
LC_COLLATE=en_CA.UTF-8
[5] LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8 LC_PAPER=en_CA.UTF-8
LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] compiler_4.1.2 fastmap_1.1.0 htmltools_0.5.2 tools_4.1.2 yaml_2.2.1
rmarkdown_2.11 knitr_1.36
[8] xfun_0.28 digest_0.6.28 rlang_0.4.12 evaluate_0.14
>
Hope this helps,
JN
On 2021-12-02 5:50 a.m., Labone, Thomas wrote:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# This works
n <- 1000# OK <= 4095
Z <- qnorm(ppoints(n))
k <- sort(rlnorm(n,log(2131),log(1.61)) / rlnorm(n,log(355),log(1.61)))
quantile(k,probs=c(0.025,0.5,0.975))
summary(k)
fit <- lm(log(k) ~ Z)
summary(fit)
gm <- exp(coef(fit)[1])
gsd <- exp(coef(fit)[2])
gm
gsd
plot(Z,k,log="y",xlim=c(-4,4),ylim=c(0.1,100))
lines(Z,gm*gsd^Z,col="red")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#this does not
n <- 10000# fails >= 4096 = 2^12
Z <- qnorm(ppoints(n))
k <- sort(rlnorm(n,log(2131),log(1.61)) / rlnorm(n,log(355),log(1.61)))
quantile(k,probs=c(0.025,0.5,0.975))
summary(k)
fit <- lm(log(k) ~ Z)
summary(fit)
gm <- exp(coef(fit)[1])
gsd <- exp(coef(fit)[2])
gm
gsd
plot(Z,k,log="y",xlim=c(-4,4),ylim=c(0.1,100))
lines(Z,gm*gsd^Z,col="red")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
