[R] Solving the equation using uniroot

[Berend Hasselman]

On 24-03-2012, at 05:35, Shant Ch wrote:

> Hello all,
> 
> I was going to solve (n-m)! * (n-k)! = 0.5 *n! * (n-m-k)!
> 
> for m when values of n and k are provided
> 
> 
> n1<-c(10,13,18,30,60,100,500)         # values of n
> 
> kx<-seq(1,7,1);                               # values of k
> 
> 
> slv2<-lapply(n1,function(n){
> 
>    slv1<-lapply(kx,function(k){
>             
> lhs<-function(m)
>              {
>                 fnk<-factorial(n-m)*factorial(n-k)-
>                       0.5*factorial(n-m-k)*factorial(n);
>                 return(fnk);
>              }
>              un2<-uniroot(lhs,c(0,n))
>              un1<-un2$root/n;
>            return(rbind(un1));
>     });
>     rjbk<-data.frame(do.call(cbind,slv1));
>     return(cbind(n,rjbk));
> });
> rj<-data.frame(do.call(rbind,slv2));
> rj
> 
> I used the above code, but I am getting this error, can anyone help me with this.
> 
> 
> Error in uniroot(lhs, c(0, n)) : f.upper = f(upper) is NA
> In addition: Warning message:
> In gamma(x + 1) : NaNs produced
> 


Don't solve   (n-m)! * (n-k)! = 0.5 *n! * (n-m-k)!
but rewrite this to use logarithms.

Use lfactorial instead of factorial to 

  fnk <- lfactorial(n-m) + lfactorial(n-k)-
              (log(0.5) + lfactorial(n-m-k) + lfactorial(n))

You will still get warnings from lfactorial:

 In lfactorial(n - m - k) : value out of range in 'lgamma'


but it's up to you to do something about that.

Berend