Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Approximating GCF
Replies: 6   Last Post: Aug 1, 2014 8:05 AM

 Messages: [ Previous | Next ]
 James Waldby Posts: 545 Registered: 1/27/11
Re: Approximating GCF
Posted: Jul 21, 2014 9:56 AM

On Mon, 21 Jul 2014 04:34:30 -0700, ryan.mitchley wrote:

> Thanks for your suggestions, James - firstly for stating my problem better than I had (I suppose that this is half of the battle), and then for the continued fraction approach. This seems like a good basis.
>
> I also had thoughts of Fourier methods (having used a cepstral algorithm for echo analysis a long time ago), but it seemed an unusual approach to take.

I think the continued fraction approach can produce an unbiased estimate
of the interval for a set of data, which is a strength; but as drawbacks
it has the issues of (a) there are O(n^4) ratios (from n choose 4) for n
data values, and (b) a clustering or measure-of-quality calculation is
needed to decide which estimates of intervals should be averaged. Issue
(a) is the reason I suggested taking a random selection of ratios to get
trial values of r -- for example, draw a few hundred or few thousand random
4-tuples from the data, and compute r for each; then use an iterative
method to gradually shift r and w so as to improve some quality measure.

BTW, here is the code that produced the array of intervals in my previous
post. Note, it has an ad hoc parameter called highjump that controls exit
from the getConvergent routine. With n and d for numerator and denominator,
if new n+d is more than highjump * old n+d, exit and return old n & d. For
example, if convergent 1/2 is followed by 26/53, return 1 and 2.

#!/usr/bin/env python
# Re 22027.msg: approximating an interval r for a data model
# yi = r * xi + w + di, where yi are observed values and other
# values are unknowns. w is an offset, di is random noise.

def getConvergent(h): # Return a continued fraction convergent of h
# highjump is a threshhold to test if most of h is accounted for
c=1; dx=d=0; m=0; nx=n=1; t=h; highjump = 20
for i in range(6):
k = int(t)
v=m; m=n; n=v+k*n
v=c; c=d; d=v+k*d
t = t - k
if (n+d) > highjump*(nx+dx) : break
dx=d; nx=n
if t < 1e-6: break
t=1/t
return (nx,dx)

def getRatio(p,q,r,s):
da = q-p
db = s-r
(da,db) = (da,db) if da < db else (db,da)
r = da/db
print '{:8.4f} {:8.4f} {:8.4f} '.format(da, db, r)
(n,d) = getConvergent(r)
if n:
return [da/n, db/d]
else:
return []

s=sorted([-1.3375, 2.0776, -4.7121, 2.0996, 3.8090])
L=len(s)
vs = []
for i in range(L):
q = s[0:i]+s[i+1:] # Get a quad
vs += getRatio(q[0],q[1],q[2],q[3])
vs += getRatio(q[0],q[2],q[1],q[3])
vs += getRatio(q[1],q[2],q[0],q[3])
print sorted([round(x,5) for x in vs])

def shoAverage(vs1, vs2):
scal = sum(vs1+vs2)/(len(vs1)+2*len(vs2))

shoAverage([x for x in vs if 1<x<2], [])
shoAverage([], [x for x in vs if 2<x])
shoAverage([x for x in vs if 1<x<2], [x for x in vs if 2<x])

--
jiw

Date Subject Author
7/17/14 James Waldby
7/20/14 James Waldby
7/21/14 Ryan M
7/21/14 James Waldby
8/1/14 Ryan M