Approximate density from moments and quantiles, then sample from it - Cross Validated most recent 30 from stats.stackexchange.com 2019-09-22T17:08:46Z https://stats.stackexchange.com/feeds/question/403451 https://creativecommons.org/licenses/by-sa/4.0/rdf https://stats.stackexchange.com/q/403451 5 Approximate density from moments and quantiles, then sample from it ilprincipe https://stats.stackexchange.com/users/9030 2019-04-16T21:01:11Z 2019-04-24T19:24:36Z <p><strong>Situation</strong></p> <p>I need to send R code to a third party to run estimations for me (I will not be able to work with the data directly). I want to simulate data to test some of the estimators before sending the code to them. </p> <p>The data provider has given me the following summary statistics for all variables: The first four moments (mean, variance, skewness, kurtosis), the four largest/smallest values (min/max), and the 1%, 5%, 10%, 25%, 50%, 75%, 90%, 95%, 99% percentiles.</p> <p>I believe the output was obtained in Stata, I am pasting an example below. The example is not from the actual data (actual data has close to 6 million observations). </p> <p><strong>Question</strong></p> <p>What would be the best way to simulate data? Initially I was just going to pick a distribution and sample from it (e.g. binomial/multinomial/(log-)normal/exponential/truncated normal), but with the information provided I presumed it is possible to do better, at least for variables which aren't binomial.</p> <p><strong>Input example</strong></p> <pre><code> Percentiles Smallest 1% 163 99 5% 216 111 10% 248 113 Obs 3,170 25% 322 114 Sum of Wgt. 3,170 50% 494 Mean 1262.359 Largest Std. Dev. 3093.165 75% 984 41584 90% 2450 54413 Variance 9567670 95% 5181 58477 Skewness 10.59025 99% 10826 59349 Kurtosis 157.7004 </code></pre> <p><strong>What I am currently doing</strong></p> <p>I have tried fitting the parameters of a specific distribution I think might be appropriate using functions in the <code>library(rriskDistributions)</code> (e.g. <code>get.lnorm.par()</code>). This works well sometimes but often it does not.</p> <p>Currently I am fitting the CDF using splines, obtain the PDF using the spline functions derivative, and then sample from it.</p> <p>Neither of these approaches work very well in general. I'm hoping for an approach that is generic and delivers a good approximation without me manually eyeballing the distribution and investigating the fit. I understand that this may be asking a lot given the limited data at my disposal.</p> <pre class="lang-r prettyprint-override"><code>## function for spline interpolation splsample &lt;- function(p, v, size = 1000000, vmin = min(v), vmax = max(v), gridsize = min(3*(vmax-vmin), 1000), step = NULL, plot = FALSE, ...) { s &lt;- splinefun(v, p, ...) if(is.null(step)) { grid &lt;- seq(from = vmin, to = vmax, length.out = gridsize) } else { grid &lt;- seq(from = vmin, to = vmax, by = step) } pr &lt;- s(grid, deriv = 1) pr[pr &lt; 0] &lt;- 0 if (plot == TRUE) { plot(grid, pr) } bs &lt;- sample(grid, p = pr, size = size, replace = TRUE) return(bs) } ## input percentiles &lt;- c(0.01, 0.05, 0.10, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99) values &lt;- c(163, 216, 248, 322, 494, 984, 2450, 5181, 10826) ## spline approximation of pdf x &lt;- splsample(percentiles, values, plot = TRUE) summary(x) mean(x) var(x) ## alternative: fitting a truncated normal library("rriskDistributions") library("msm") dpar &lt;- get.tnorm.par(p = percentiles, q = values) x &lt;- rtnorm(10000, mean = dpar["mean"], sd = dpar["sd"], lower = dpar["lower"], upper = dpar["upper"]) x[x &lt; 0] &lt;- 0 summary(x) mean(x) var(x) </code></pre> https://stats.stackexchange.com/questions/403451/-/404882#404882 1 Answer by Anthony for Approximate density from moments and quantiles, then sample from it Anthony https://stats.stackexchange.com/users/27148 2019-04-24T19:19:09Z 2019-04-24T19:24:36Z <p>To quickly simulate based on moments, try <code>rpearson()</code> from <code>library(PearsonDS)</code>.</p> <pre class="lang-r prettyprint-override"><code>library(PearsonDS) target.moms &lt;- c(1262.39, 9567670, 10.59025, 157.7004) y &lt;- rpearson(n=1000000, moments=target.moms) </code></pre> <p><code>rpearson()</code> works well for matching the moments. However, the splines approach that you're already using will be better at recovering the percentiles. See below for an example.</p> <pre class="lang-r prettyprint-override"><code>#Evaluating the results library(moments) eval &lt;- function(data) { result.list &lt;- list(mean=mean(data), var=var(data), skew=skewness(data), kurt=kurtosis(data), quantile(data, c(.01,.05,.10,.25,.50,.75,.90,.95,.99) ) ) round(unlist(result.list), 2) } x &lt;- splsample(percentiles, values, plot = TRUE) eval(x) #splines eval(y) #rpearson() </code></pre>