[R] Help in R

[Charles C. Berry]

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 <djmuser at 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.hlx at 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 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<http://www.r-project.org/posting-guide.html>
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>
>>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> 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.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901