Solving Systems of Non-linear equations
Another follow up comment. I tried it in Maxima (also free) and noticed that it has the capability of performing the solution in just a single line using the Maxima solve command so you may prefer that. Note that the first line display2d:false turns off fancy 2d output and you can omit it if you want the fancy 2d output. display2d:false; solve([mean = a/(a+b), variance = (a*b)/(((a+b)^2) * (a+b+1))], [a,b]); The output looks like this: (%i18) display2d:false; (%o18) FALSE (%i19) solve([mean = a/(a+b), variance = (a*b)/(((a+b)^2) * (a+b+1))], [a,b]); (%o19) [[a = -(mean*variance+mean^3-mean^2)/variance, b = ((mean-1)*variance+mean^3-2*mean^2+mean)/variance],[a = 0,b = 0]]
On 11/30/05, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
Just one addition to this. I noticed that its not really true that the output can be used in R verbatim since the C output uses pow instead of ^; however, if one replaces the "code c" statement with the statement "list export" then it is valid R. That is the input to mathomatic should be: mean = a/(a+b) variance = (a*b)/(((a+b)^2) * (a+b+1)) eliminate b a simplify list export eliminate a b simplify list export On 11/30/05, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
Sorry I seemed to have messed up the copying and pasting. Here it is again. --- Go to http://mathomatic.orgserve.de/math/ and install mathomatic (its free) or just connect to the online server and do this. The C output, i.e the result of the two code c commands, can be used verbatim in R. Note that mathomatic does not support logs but for simple problems like this its very useful. Note that 1-> and 2-> are the mathomatic prompts and what comes after them are what I typed in. The entry goes into the corresponding equation space, i.e. equation 1 or equation 2. This is what you enter: mean = a/(a+b) variance = (a*b)/(((a+b)^2) * (a+b+1)) eliminate b a simplify code c eliminate a b simplify code c and this is the entire session: 1-> mean = a/(a+b) a #1: mean = ------- (a + b) 1-> variance = (a*b)/(((a+b)^2) * (a+b+1)) a*b #2: variance = ------------------------- (((a + b)^2)*(a + b + 1)) 2-> eliminate b Solving equation #1 for (b)... 1 (a^2)*(---- - 1) mean #2: variance = --------------------------------------------------- 1 1 (((a + (a*(---- - 1)))^2)*(a + (a*(---- - 1)) + 1)) mean mean 2-> a mean*(1 - mean) #2: a = mean*(--------------- - 1) variance 2-> simplify ((mean^2) - (mean^3)) #2: a = --------------------- - mean variance 2-> code c a = ((((mean * mean) - pow(mean, 3.0)) / variance) - mean); 2-> eliminate a Solving equation #1 for (a)... b*mean ((mean^2) - (mean^3)) #2: ---------- = --------------------- - mean (1 - mean) variance 2-> b mean*(1 - mean) #2: b = (--------------- - 1)*(1 - mean) variance 2-> simplify ((mean^2) - mean) #2: b = (1 + -----------------)*(mean - 1) variance 2-> code c b = ((1.0 + (((mean * mean) - mean) / variance)) * (mean - 1.0)); On 11/30/05, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
Go to http://mathomatic.orgserve.de/math/ and install mathomatic (its free) or just connect to the online server and do this. The C output, i.e the result of the two code c commands, can be used verbatim in R. Note that mathomatic does not support logs but for simply problems like this its very useful. Note that 1-> and 2-> are the mathomatic prompts and what comes after them are what I typed in. The entry goes into the corresponding equation space, i.e. equation 1 or equation 2. 1-> mean = a/(a+b) a #1: mean = ------- (a + b) 1-> variance = (a*b)/(((a+b)^2) * (a+b+1)) a*b #2: variance = ------------------------- (((a + b)^2)*(a + b + 1)) 2-> eliminate b Solving equation #1 for (b)... 1 (a^2)*(---- - 1) mean #2: variance = --------------------------------------------------- 1 1 (((a + (a*(---- - 1)))^2)*(a + (a*(---- - 1)) + 1)) mean mean 2-> a mean*(1 - mean) #2: a = mean*(--------------- - 1) variance 2-> simplify ((mean^2) - (mean^3)) #2: a = --------------------- - mean variance 2-> eliminate a Solving equation #1 for (a)... b*mean ((mean^2) - (mean^3)) #2: ---------- = --------------------- - mean (1 - mean) variance 2-> b mean*(1 - mean) #2: b = (--------------- - 1)*(1 - mean) variance 2-> simplify ((mean^2) - mean) #2: b = (1 + -----------------)*(mean - 1) variance 2-> code c b = ((1.0 + (((mean * mean) - mean) / variance)) * (mean - 1.0)); 2-> #1 b*mean #1: a = ---------- (1 - mean) 1-> code c a = (b * mean / (1.0 - mean)); On 11/30/05, Scott Story <sstory at montana.edu> wrote:
I am trying to write a function that will solve a simple system of nonlinear equations for the parameters that describe the beta distribution (a,b) given the mean and variance. mean = a/(a+b) variance = (a*b)/(((a+b)^2) * (a+b+1)) Any help as to where to start would be welcome. -- Scott Story Graduate Student MSU Ecology Department 319 Lewis Hall Bozeman, Mt 59717 406.994.2670 sstory at montana.edu
______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html