Hi Marc
I think that we could help you better if we knew in which context you need
sample from a sum constrained normal distribution. However this is more a
question on probability theory than on how to do it in R.
The proposal so far has been linear transformation of multivariate normal
distribution (Marc, Rui), mixture of normal and reflected normal distribution
(Boris, try that with e.g. mu = 2), normal distribution mixed with single point
with positive mass (Jlucke), degenerated normal distribution (Greg).
What you in fact want to do is to draw samples from a conditional distribution.
The condition is the sum constraint so if we have x = (x1, x2, ..., xn) then
sum_{i=1}^n xi = 0 or x1 + x2 + ... x{n-1} = xn so you want to draw samples
from P(x given that x is normal distributed and sum(x)=0). The sum constraint
gives in fact what is called distributions on the simplex. Google for "normal
distribution simplex" and you will get almost 2 mill hits. The second shows how
to sample using Gibbs sampling
(http://dobigeon.perso.enseeiht.fr/papers/Dobigeon_TechReport_2007b.pdf).
However you can probably just use other distributions given sum constraint
since you say that you only need the sample as initial values for a MCMC
algorithm. Many methods are available from "compositional statistics" (google
for that, Aitchison 1986 is the pioneer).
At least two packages are available for R:"compositions" with the latest
version from 2013 and can be found in the archives and "robComposition" still
maintained.
Hope that helps, it is help for yourself to find a solution.
Yours sincerely / Med venlig hilsen
Frede Aakmann Tøgersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling
Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
fr...@vestas.com
http://www.vestas.com
Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender.
> -----Original Message-----
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> On Behalf Of Marc Marí Dell'Olmo
> Sent: 1. april 2014 16:57
> To: Boris Steipe
> Cc: r-help@r-project.org
> Subject: Re: [R] A vector of normal distributed values with a sum-to-zero
> constraint
>
> Boris is right. I need this vector to include as initial values of a
> MCMC process (with openbugs) and If I use this last approach sum(x)
> could be a large (or extreme) value and can cause problems.
>
> The other approach x <- c(x, -x) has the problem that only vectors
> with even values are obtained.
>
> Thank you!
>
>
> 2014-04-01 16:25 GMT+02:00 Boris Steipe <boris.ste...@utoronto.ca>:
> > But the result is not Normal. Consider:
> >
> > set.seed(112358)
> > N <- 100
> > x <- rnorm(N-1)
> > sum(x)
> >
> > [1] 1.759446 !!!
> >
> > i.e. you have an outlier at 1.7 sigma, and for larger N...
> >
> > set.seed(112358)
> > N <- 10000
> > x <- rnorm(N-1)
> > sum(x)
> > [1] -91.19731
> >
> > B.
> >
> >
> > On 2014-04-01, at 10:14 AM, jlu...@ria.buffalo.edu wrote:
> >
> >> The sum-to-zero constraint imposes a loss of one degree of freedom. Of
> N samples, only (N-1) can be random. Thus the solution is
> >> > N <- 100
> >> > x <- rnorm(N-1)
> >> > x <- c(x, -sum(x))
> >> > sum(x)
> >> [1] -7.199102e-17
> >>
> >> >
> >>
> >>
> >>
> >>
> >>
> >>
> >>
> >>
> >> Boris Steipe <boris.ste...@utoronto.ca>
> >> Sent by: r-help-boun...@r-project.org
> >> 04/01/2014 09:29 AM
> >>
> >> To
> >> Marc Marí Dell'Olmo <marceivi...@gmail.com>,
> >> cc
> >> "r-help@r-project.org" <r-help@r-project.org>
> >> Subject
> >> Re: [R] A vector of normal distributed values with a sum-to-zero
> constraint
> >>
> >>
> >>
> >>
> >>
> >> Make a copy with opposite sign. This is Normal, symmetric, but no longer
> random.
> >>
> >> set.seed(112358)
> >> x <- rnorm(5000, 0, 0.5)
> >> x <- c(x, -x)
> >> sum(x)
> >> hist(x)
> >>
> >> B.
> >>
> >> On 2014-04-01, at 8:56 AM, Marc Marí Dell'Olmo wrote:
> >>
> >> > Dear all,
> >> >
> >> > Anyone knows how to generate a vector of Normal distributed values
> >> > (for example N(0,0.5)), but with a sum-to-zero constraint??
> >> >
> >> > The sum would be exactly zero, without decimals.
> >> >
> >> > I made some attempts:
> >> >
> >> >> l <- 1000000
> >> >> aux <- rnorm(l,0,0.5)
> >> >> s <- sum(aux)/l
> >> >> aux2 <- aux-s
> >> >> sum(aux2)
> >> > [1] -0.000000000006131392
> >> >>
> >> >> aux[1]<- -sum(aux[2:l])
> >> >> sum(aux)
> >> > [1] -0.00000000000003530422
> >> >
> >> >
> >> > but the sum is not exactly zero and not all parameters are N(0,0.5)
> >> > distributed...
> >> >
> >> > Perhaps is obvious but I can't find the way to do it..
> >> >
> >> > Thank you very much!
> >> >
> >> > Marc
> >> >
> >> > ______________________________________________
> >> > R-help@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.
> >>
> >> ______________________________________________
> >> R-help@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.
> >>
> >
> > ______________________________________________
> > R-help@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.
>
> ______________________________________________
> R-help@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.
______________________________________________
R-help@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.
