Re: [R] Mixture of Distributions

From: Ted Harding <Ted.Harding_at_manchester.ac.uk>
Date: Thu, 21 Feb 2008 22:20:07 +0000 (GMT)


I just realised I made a bad mistaqke (see below)

On 21-Feb-08 21:39:56, Ted Harding wrote:
> On 21-Feb-08 20:58:25, Evgenia wrote:

>> 
>> Dear R users, 
>> I would like to sample from a mixture distribution
>> p1*f1+p2*f2+p3*f3 with f1,f2,f3 three different forms
>> of distributions. I know that in the case of two
>> distributions I have to sample the mixture compoment
>> membership.
>> 
>> How can I weight my three distributions with their
>> respective probabilities?
>> 
>> Evgenia

>
> There are several ways in which you could write code
> to do this, but all amount to the fact that each value
> sampled from such a mixture is obtained by
> a) choose between f1, f2, f3 at random according to the
> probabilities p1, p2, p3 (it is assumed p1+p2+p3=1).
> b) sample 1 value from whichever of f1, f2, f3 was chosen.
>
> Suggestion:
>
> Suppose the functions rf1(n), rf2(n), rf3(n) respectively
> sample n values from f1, f2 and f3.

S0 <- cbind(rf1(n),rf2(n),rf3(n))
ix <- sample(c(1,2,3),n,3,prob=c(p1,p2,p3),replace=TRUE)

So far so good!

S <- S0[,ix]

But this will not do what I intended (i.e. select element ix[1] from the first row os S0, ix[2] from the second row, and so on).

Instead, it will return an nxn matrix with n rows, and n columns which are copies of columns of S0 in the order selected by ix!

The following will do the correct thing, though there must be a neater way of doing it! The example below has been corrected.

S <- S0[,1]*(ix==1) + S0[,2]*(ix==2) + S0[,3]*(ix==3)

> will produce n values, each of which is a sample of size 1
> from the mixture.
>
> Example:

rf1 <- function(n){rnorm(n,0,1)} ## normal distribution with mean=0
rf2 <- function(n){rnorm(n,4,1)} ## normal distribution with mean=4
rf3 <- function(n){rnorm(n,9,1)} ## normal distribution with mean=9

p1 <- 0.2; p2 <- 0.3; p3 <-0.5

S0 <- cbind(rf1(500),rf2(500),rf3(500))
ix <- sample(c(1,2,3),500,prob=c(p1,p2,p3),replace=TRUE) S <- S0[,1]*(ix==1) + S0[,2]*(ix==2) + S0[,3]*(ix==3)

hist(S,breaks=50)

>
> Hoping this helps.
> Ted.

And again!
Ted.



E-Mail: (Ted Harding) <Ted.Harding_at_manchester.ac.uk> Fax-to-email: +44 (0)870 094 0861
Date: 21-Feb-08                                       Time: 22:20:04
------------------------------ XFMail ------------------------------

______________________________________________
R-help_at_r-project.org mailing list
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. Received on Thu 21 Feb 2008 - 22:23:32 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 04 Mar 2008 - 04:30:18 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive