
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 measureofquality 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 4tuples 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 < 1e6: break t=1/t return (nx,dx)
def getRatio(p,q,r,s): da = qp db = sr (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

