[R] Solving Equations

[Gabor Grothendieck]

On Sun, Jan 22, 2012 at 1:25 PM, Eliano <eliano.m.marques at gmail.com> wrote:
> People,
>
> I'm researching some Bayesian statistic topics and in the midle of my study
> i found a very simple problem and i'm trying to find a simple package to
> solve this type of equations:
>
> Lets say that i need to compute beta values for the beta distribution and i
> now for example:
>
> E(teta)=a/(a+b) = 0,5
> Var(teta)=ab/((a+b)^2(a+b+1))=0.05
>
> So if i want to solve this to non-linear system to find a,b for the beta
> distribution wich pack should i use?
>

The equations giving a and b as a function of m and v are:

> library(Ryacas)
> a <- Sym("a"); b <- Sym("b")
> m <- Sym("m"); v <- Sym("v")
> Solve( List(a/(a+b) == m, a*b/((a+b)^2*(a+b+1)) == v), List(a, b) )
expression(list(list(a == m^2 * (1 - m)/v - m, b == a/m - a)))

Based on the above we write this R function:

beta.parms <- function(m, v) {
   a <- m^2 * (1-m)/v - m
   b <- a/m - a
   c(a = a, b = b)
}

and run it:

> beta.parms(m = 0.5, v = 0.05)
a b
2 2

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com