Thanks for the hint. I wrote the following:
> n <- 10
> a <- 1
> b<- 2
> set.seed(1)
> y <- rnorm(n)
>
> xmat <- as.matrix( expand.grid( rep( list(0:1), n ) ))
>
> s <- numeric(2^(n))
> for (i in 1:2^(n)){
+ s[i]<- beta(a+rowSums(xmat)[i], b+n-rowSums(xmat)[i])*
+ prod((1-xmat[i,])*dnorm(y,0,1)+xmat[i,]*dnorm(y,0,1))}
> sum(s)
[1] 3.015227e-06
There is still a problem though. when n is large, "expand.grid" seems not to
work.
> expand.grid( rep( list(0:1), 100 ) )
Error in rep.int(rep.int(seq_len(nx), rep.int(rep.fac, nx)), orep) :
invalid 'times' value
In addition: Warning message:
In rep.int(rep.int(seq_len(nx), rep.int(rep.fac, nx)), orep) :
NAs introduced by coercion
Is there a way to correct it?
Thank you.
2009/12/15 Charles C. Berry <cbe...@tajo.ucsd.edu>
> On Tue, 15 Dec 2009, li li wrote:
>
> Hi:
>> Thank you for your reply. I realize that I did not state the problem
>> clearly before.
>> Here is the problem again.
>>
>> Let X = (X1, X2, ..., Xn) and Y=(Y1,Y2,...,Yn). Xi's can be 0 or 1.
>> When Xi=1, Yi is distributed as N(2,1). When Xi=0, Yi is distributed as
>> N(0,1).
>> There are 2^n possible X values. For example, x=(0,0,...,0) , i.e. all xi
>> is
>> 0.
>>
>> The summand is:
>>
>> beta(a+sum(xi), b+n-sum(xi) ) * [ (1-x1) * dnorm(y1, 0, 1) + x1 *
>> dnorm(y1,
>> 2, 1) ] *
>> [(1-x2) * dnorm(y2, 0, 1) + x2 * dnorm(y2, 2, 1) ] * ...* [ (1-xn) *
>> dnorm(yn, 0, 1) + xn * dnorm(yn, 2, 1) ]
>>
>> where sum(xi)=x1 + x2 + ... + xn.
>>
>> For example, when x=(0,0,...,0) , i.e. all xi is 0. The summand is
>>
>> beta(a, b+n) * dnorm(y1, 0, 1) * dnorm(y2, 0, 1) * ...* dnorm(yn, 0, 1)
>>
>>
>> I want to take the sum of the summand over all possible X values. There
>> are
>> 2^n of them.
>>
>> I am wondering whether there is a function in R that can take care of this
>> kind of sum.
>>
>>
> A function that does exactly that? I do not know of one, but have not
> searched very hard.
>
> If I understood what you want (and "commented, minimal, self-contained,
> reproducible code" was requested and would have helped me), a simple one
> line command would give the result.
>
> For
>
> set.seed(1)
>> n <- 10
>> y <- rnorm(n)
>> x <- rbinom(10,1,.5)
>>
>> a <- 1
>> b <- 2
>>
>
>
> The answer is given by:
>
>> # [ one line command omitted ]
>>
> [1] 1.298587e-08
>
>>
>>
> You will need to understand "vectorization" and "recycling rules" to do
> this neatly in one line.
>
> It will help you to study the sections of the manuals that pertain to those
> topics (and perhaps the mail ist archive) as well as to study the help and
> examples for
>
> ?dnorm
> ?beta
> ?expand.grid
> ?matrix
> ?rowSums
> ?as.matrix
>
>
> And here is a helpful hint
>
> xmat <- as.matrix( expand.grid( rep( list(0:1), n ) )
>
> will be useful.
>
> If you need more help, be sure to review the Posting Guide and run the
> function help.request() to help you better frame your question.
>
> HTH,
>
> Chuck
>
> p.s. Is this homework? If so, ask your instructor for help first.
>
>
> Thank you
>>
>>
>>
>> 2009/12/15 Dennis Murphy <djmu...@gmail.com>
>>
>> Hi:
>>>
>>> Your expression makes no sense in that I don't see where the summation
>>> can
>>> obtain 'over all 2^n possible values...' the way you wrote it. How about
>>> this?
>>>
>>> Let x = (x1, x2, ..., xn). The summand can be then be expressed in R as
>>>
>>> beta(a + sum(x), b + n - sum(x)) * ((1 - x) * dnorm(x, 0, 1) + x *
>>> dnorm(x, 2, 1)) .
>>>
>>> The beta function term is a scalar, the second expression is a vector.
>>> Do you want the sum of the products of the vector elements? If so, it
>>> should be as easy as
>>>
>>> beta(a + sum(x), b + n - sum(x)) * sum((1 - x) * dnorm(x, 0, 1) + x *
>>> dnorm(x, 2, 1))
>>>
>>> If, on the other hand, you meant to write the summand as
>>>
>>> beta(a + x, b + 1 - x) * ((1 - x) * dnorm(x, 0, 1) + x * dnorm(x, 2, 1))
>>> ,
>>>
>>> then the sum is
>>>
>>> sum(beta(a + x, b + 1 - x) * ((1 - x) * dnorm(x, 0, 1) + x * dnorm(x, 2,
>>> 1))
>>>
>>>
>>> More below...HTH
>>> Dennis
>>>
>>> On Mon, Dec 14, 2009 at 6:38 PM, li li <hannah....@gmail.com> wrote:
>>>
>>> Hello,
>>>> Can anyone give me some suggestion in term of calculating the sum
>>>> below.
>>>> Is there a function in R that can help doing it faster?
>>>>
>>>> x1, x2, ...xn where xi can be 0 or 1. I want to calculate the following:
>>>>
>>>> sum{ beta[a+sum(xi), b+n-sum(xi) ]* [ (1-x1)dnorm(0,1)+x1dnorm(2,1) ]*
>>>> [
>>>> (1-x2)dnorm(0,1)+x2dnorm(2,1) ]* ...* [ (1-xn)dnorm(0,1)+xndnorm(2,1) ]
>>>> }
>>>>
>>>>
>>> The problem I have is that you're taking the product of the n normal
>>> mixtures and then
>>> you somehow want to sum over them; if
>>> this is meant to be a likelihood function, then it seems that you would
>>> want
>>>
>>> beta(a+sum(x), b+n-sum(x)) * prod((1 - x) * dnorm(x, 0, 1) + x *
>>> dnorm(x,
>>> 2, 1))
>>>
>>> instead. As I said up front, it's not clear what you really want, so I
>>> leave you with
>>> multiple possibilities in the hope that one of them is what you need...
>>>
>>>
>>> The sum in the beginning is over all 2^n possible values for the
>>>> vector
>>>> x1,
>>>> x2, ...xn .
>>>>
>>>> Thank you very much!
>>>>
>>>> Hannah
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> ______________________________________________
>>>> 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<http://www.r-project.org/posting-guide.html>
>>>> <http://www.r-project.org/posting-guide.html>
>>>>
>>>> and provide commented, minimal, self-contained, reproducible code.
>>>>
>>>>
>>>
>>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> 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<http://www.r-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
> Charles C. Berry (858) 534-2098
> Dept of Family/Preventive
> Medicine
> E mailto:cbe...@tajo.ucsd.edu UC San Diego
> http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
>
>
>
[[alternative HTML version deleted]]
______________________________________________
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.
