Re: [R] Solving Systems of Non-linear equations

From: Gabor Grothendieck <ggrothendieck_at_gmail.com>
Date: Tue 06 Dec 2005 - 02:31:12 EST

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@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@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@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@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@montana.edu
> > > >
> > > > ______________________________________________
> > > > R-help@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
> > > >
> > >
> >
>



R-help@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 Received on Tue Dec 06 02:56:25 2005

This archive was generated by hypermail 2.1.8 : Tue 06 Dec 2005 - 04:25:48 EST