[R] solving equation consisting of factorials

David Winsemius dwinsemius at comcast.net
Sat Oct 30 18:09:06 CEST 2010


On Oct 30, 2010, at 9:55 AM, Shant Ch wrote:

>
>
> Hi,
>
>
> I was trying to solve for m such that,
>
> choose(n-1,k-1)+choose(n-2,k-1)+...+choose(n-m+1,k-1) - ck=0
>
> ck is calculated by the way given in the code. Now I solved for m  
> using uniroot
> but the error messages appeared.
>
> n<-30
> k=3
> alpha<-0.05
> nk<-choose(n,k);
> ck<-(nk/2)+k*nk*qnorm(1-alpha)/(12*n)    # computation of ck.

Given this obviously non-integer value for ck and the obviously  
integer sum of the preceding terms in your equation above, why should  
we think there would be ANY solution?


>
> lhs<-function(z)
> {
>    Smc=0;
>    fnk=0
>    for(y in 1:z)
>    {
>       fnk[y]<-choose((n-y+1),(k-1))
>       Smc=Smc+fnk[y]-ck
>    }
>    return(Smc);
> }
>
>
> In order to solve the equation first I tried to find out the point  
> near which
> lhs=0 using curve function.
>
> curve(lhs,-70,-60)

Try instead:

 > plot(-60:70, sapply(-60:70, lhs) )

Which will show you that there are no values of lhs() in the range you  
are exploring , probably because you are trying to index that y object  
with negative values inside that lhs function.

Maybe you should describe what your real goal of all this is. There  
are packages for integer programing, but you need to provide a better  
problem specification ( .i.e, a problem that is solvable.)


>
> Error in xy.coords(x, y, xlabel, ylabel, log) :
>  'x' and 'y' lengths differ
> In addition: Warning message:
> In 1:z : numerical expression has 101 elements: only the first used
>
>
> Then I computed lhs function for every z  one by one, I found that  
> lhs changes
> its sign between (-64,-63), so I used the following uniroot function.

Must have been a different function:
 > lhs(-63)
numeric(0)

-- 
David.
>
> uniroot(lhs(-64,-63))
>
> Error in min(interval) : 'interval' is missing
>
>
> Can anyone help me out? Thanks.
>
> Shant


David Winsemius, MD
West Hartford, CT



More information about the R-help mailing list