# Re: [R] Fit model to data and use model for data generation

From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>
Date: Thu 25 Jan 2007 - 08:25:11 GMT

On Wed, 24 Jan 2007, Stephen D. Weigand wrote:

>
> On Jan 24, 2007, at 10:34 AM, Benjamin Otto wrote:
>
>> Hi,
>>
>> Suppose I have a set of values x and I want to calculate the
>> distribution of
>> the data. Ususally I would use the "density" command. Now, can I use
>> the
>> resulting "density-object" model to generate a number of new values
>> which
>> have the same distribution? Or do I have to use some different
>> function?
>>
>> Regards,
>>
>> Benjamin
>>
>> --
>> Benjamin Otto
>> Universitaetsklinikum Eppendorf Hamburg
>> Institut fuer Klinische Chemie
>> Martinistrasse 52
>> 20246 Hamburg
>>
>
> You could sample from the x's in the density object with probability
> given by the y's:

That gives a discrete distribution, which may well matter for small samples.

x.new <- rnorm(n, sample(x, size = n, replace=TRUE), bw)

where bw is the bandwidth used by density (d\$bw in this example). (This is known as a 'smoothed bootstrap' in some circles.)

> ### Create a bimodal distribution
> x <- c(rnorm(25, -2, 1), rnorm(50, 3, 2))
> d <- density(x, n = 1000)
> plot(d)
>
> ### Sample from the distribution and show the two
> ### distributions are the same
> x.new <- sample(d\$x, size = 100000, # large n for proof of concept
> replace = TRUE, prob = d\$y/sum(d\$y))
> dx.new <- density(x.new)
> lines(dx.new\$x, dx.new\$y, col = "blue")

BTW, lines(density(x.news), col = "blue") works here, and you do need to remember that a kde is biased. But my solution matches better than yours.

```--
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help