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
Scott Story <sstory at montana.edu> writes:> 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.On a pad of paper... First look at mean*(1-mean)/variance, and the rest should follow. (Hint minimized in case this was homework...) -- O__ ---- Peter Dalgaard ??ster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
On 11/30/2005 10:14 AM, Scott Story 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.You should use a package like Maple or Mathematica (or just some pencil and paper work) to determine the solution. Then the function is really easy to write. Maple gives > solve({mean = a/(a+b),variance = (a*b)/(((a+b)^2) * (a+b+1))},{a,b}); {a = -mean*(variance+mean^2-mean)/variance, b = (variance+mean^2-mean)*(mean-1)/variance} from which you can write your own function pretty easily. Duncan Murdoch
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
>
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
> >
>
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 > > > > > >
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. > >The simplest way to go is probably to use the fitdistr function in MASS package. YNOEL
>>>>> "YNoel" == NOEL Yvonnick <yvonnick.noel at uhb.fr> >>>>> on Thu, 01 Dec 2005 14:42:44 +0100 writes:YNoel> 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. >> >> YNoel> The simplest way to go is probably to use the fitdistr function in MASS YNoel> package. yes, definitely for estimating the beta parameters - namely do something smarter than the above moment estimates. However, I (and probably most respondents to this thread) believe that the question was more general (as the subject indicates) about solving non-lin eq systems, rather than about estimating the beta parameters (by a sub-optimal method). Martin Maechler, ETH Zurich
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 > > > > > > > > > >