[VA] SRC#005- April, 1st Mean Minichallenge

[Mike (Austria)]

Thanks Valentin for this inspiring minichallenge. I particularly enjoy math riddles when they involve little-known facts about well-known mathematics. This one definitely fits this description particularly well - the fact that the various means are discrete members of a continuous family had escaped my attention so far (for decades - embarrassingly enough :-) ).

Since there have been no solutions for RPL, let me share my take on this. It should work for any version of RPL starting with the HP48G(X).

1. RPL function GM for the Generalized Mean function \(G_p(x_1,\ldots,x_n)=(\frac{1}{n}\sum_{i=1}^{n}x_i^p)^{1/p}\).

Code:

«
  @  {list} real  →  real
  @      {x_i} p  →  gm(x,p)

  OVER SIZE → x p n     @ n is size(x)
  «
    x                   @ prepare x on stack
    IF p 0 == THEN      @ geometric mean
      ΠLIST             @ prod(x_i)
      n                 @ prepare exponent on stack
    ELSE                @ general case
      p ^ ΣLIST n /     @ sum(x_i^p) / n
      p
    END
    XROOT               @ % ^(1/p) or % ^(1/n)
  »
»
'GM' STO

Size: 106.5 bytes.
Exponent \(p = k-2\)
Usage: put list {\(x_i\)} and exponent \(p\) on stack, invoke GM.
I like the terseness of this function.

2. Test: four standard means

Code:

{ 4 1 20.19 } 'X' STO       @ save list
X -1 GM ->  2.30852787042   @ harmonic
X  0 GM ->  4.32247115133   @ geometric
X  1 GM ->  8.39666666667   @ arithmetic
X  2 GM -> 11.8972840038    @ quadratic

3. Solving for \(k=p+2\) such that \(G_p(x_i)=\pi\)

Auxiliary function G: for list {\(x_i\)} and exponent \(p\) in global variables X and P, return GM(x,p)-\(\pi\)

Rationale: the HP48 solver expects all parameters and the unknown to be solved for (\(p\) in this case) in global variables.

Code:

« X P GM π - »
'G' STO

Find zero of G

Code:

'G' RCL   @ recall G to stack
'P'       @ variable to solve for
0         @ initial guess
ROOT      @ solve
2 +       @ p -> k
-> 1.55425243274

4. Bonus: Code for Generalized f-Mean \(G_f(x)=f^{-1}\left(\frac{1}{n}\sum_{i=1}^{n}f(x_i)\right)\)

A further step in generalizing the mean is by replacing the \(p\)-th power of \(x_i\) by applying an arbitrary continuous strictly monotonic function and taking its inverse after computing the sum and dividing by \(n\). Not very surprising, the geometric mean results from using \(f(x)={\rm ln}(x)\) the natural logarithm. Btw, the requirement that all \(x_i\) be positive is actually not a restriction because the geometric mean is defined for positive arguments only in the first place.

Since in general the inverse of \(f\) cannot be assumed to be known analytically, we first define an auxiliary function FINV, which, given \(f, y\), and an estimate \(e\), will return an approximation of \(x=f^{-1}(y)\) by numerically solving \(f(x)=y\) for the unknown \(x\).

Code:

«
  @ function real real → real
  @ f y estimate → x such that f(x)=y

  → f y e
  «
    « tX f EVAL y - »   @ find zero of this function
    'tX'                @ temporary global variable to solve for
    e                   @ initial estimate
    ROOT
    { tX } PURGE        @ eliminate tX
  »
»
'FINV' STO

This is essentially the inverse function solver on page 2-54 of the HP48G Advanced User's Manual (AUM). I had initially hoped to use a compiled local variable (this is a local variable whose name starts with a backarrow ←) because the AUM promises that this can be used whenever a called function expects a global variable - but ROOT does not let me do this, thus the clumsy { tX } PURGE at the end.

With this complexity hidden away, it becomes straightforward to implement the generalized \(f\)-Mean function GFM:

Code:

«
  @  {list} function  →  real
  @          {x_i} f  →  gfm(x,f)


  OVER SIZE → x f n     @ n is size(x)
  «
    f                   @ prepare f: will be used by FINV
    x f EVAL            @ f(x)
    ΣLIST n /           @ sum(f(x_i)) / n
    x 1 GET             @ first element as initial guess
                        @ now on stack: f sum(%)/n e
    FINV                @ f^(-1)
  »
»
'GFM' STO

Why pick one of the original list elements to estimate the solution: we make no assumptions about the domain of \(f\). As long as the user supplies a list of valid elements, any of them is suitable as the starting point for the ROOT solver of the HP48, which is quite robust. This is cheaper than e.g. taking the arithmetic mean as an estimate or bracketing the solution by searching the minimum and maximum of {\(x_i\)}.

Examples:

Code:

X « INV » GFM -> 2.30852787043
X « LN » GFM  -> 4.32247115132
X « » GFM     -> 8.39666666667   @ empty program = identity function
X « SQ » GFM  -> 11.8972840038
X « EXP » GFM -> 19.091387809    @ smooth maximum approximation

Up to the occasional 1ULP numerical error, above results are reproduced. The last example is closely related to the LogSumExp function, which is one of the Smooth Maximum Approximations used in optimization problems.

I hope all this makes sense. Setting up these functions was fun - thanks for giving me an excuse to learn and play with stuff i'd otherwise never have looked at again. Despite RPL's deficiencies (such as handling many trivial cases incorrectly), I like the expressive power and abstraction / object orientation of this language.

P.S. This is my first post in the new forum. While I stumbled upon a nice description how to include math, I could not find a way to use a [pre] tag. So I had to reluctantly adorn my programs with code tags lest all formatting be lost. I'd gladly replace these by a better solution.