In [12]:

    

from numpy import *
from scipy import *
from sklearn import *


    

In [13]:

    

T0 = 10000
T1 = 30000


    

In [141]:

    

N = 50
alpha = 0.1
maxSVD = 0.9


    

In [44]:

    

def legendre(n, x):
    if n==0:
        return 1
    elif n==1:
        return x * (2*n+1)/2
    elif n==2:
        return (1.5*x*x-0.5) * (2*n+1)/2
    elif n==3:
        return (2.5*x*x*x-1.5*x) * (2*n+1)/2


    

In [207]:

    

a = loadtxt("input.txt")
#b2 = loadtxt("output_lin.txt")
b1 = loadtxt("output_ml_0.txt")
b2 = loadtxt("output_ml_1.txt")


    

In [208]:

    

figure()
clf()
plot(b1[1000:1400,0])
plot(b1[1000:1400,1])
plot(b1[1000:1400,2])
plot(b1[1000:1400,3])
#ylim(-1.0, 1.0)


    

    
Out[208]:





[<matplotlib.lines.Line2D at 0x10fbb6ec>]

In [209]:

    

def calc_cap(ds):
    y = zeros(T1, dtype=float)
    for i in range(T1):
        tmp = 1.0
        for j in range(len(ds)):
            tmp *= legendre(ds[j], a[T0+i-j])
        y[i] = tmp
    lin = linear_model.LinearRegression()
    lin.fit(b2[T0:T0+T1,:], y[:T1])
    z = lin.predict(b2[T0:,:])
    
    mse = 0
    yy = 0
    for i in range(T1):
        mse += (y[i] - z[i])**2
        yy += y[i]**2
    
    return 1.0 - mse/yy


    

In [210]:

    

memo = zeros(50, dtype=float)
for i in range(50):
    ds = zeros(i+1)
    ds[i] = 1
    memo[i] = calc_cap(ds)
    print i, memo[i]


    

    




0 0.975996301402
1 0.770868376788
2 0.48597818151
3 0.482204110776
4 0.385629776775
5 0.372012676998
6 0.346719779172
7 0.305546667784
8 0.297580286813
9 0.28791345431
10 0.262937112824
11 0.248051530355
12 0.246285549924
13 0.240551949135
14 0.226329833459
15 0.215795906574
16 0.211795269962
17 0.209900495986
18 0.202054021994
19 0.190713003238
20 0.182510896819
21 0.179144944824
22 0.178232464221
23 0.173741781264
24 0.1692437532
25 0.162190453111
26 0.157231961931
27 0.154597268741
28 0.153392029931
29 0.151261629521
30 0.147103576997
31 0.143248289621
32 0.140185752889
33 0.137448043465
34 0.135375472762
35 0.134403152423
36 0.132889022747
37 0.13058048352
38 0.127668087164
39 0.123931636721
40 0.120890269542
41 0.119473520583
42 0.118054084315
43 0.117259572607
44 0.115827892408
45 0.114421023675
46 0.112703430108
47 0.110758805048
48 0.108127580781
49 0.106287407489

In [211]:

    

#clf()
figure(10)
plot(memo)


    

    
Out[211]:





[<matplotlib.lines.Line2D at 0x10fe566c>]

In [19]:

    

memo2 = zeros(1000, dtype=float)
cnt = 0
for i in range(1, 10):}
    for j in range(i, 10):
        ds = zeros(10)
        ds[i] += 1
        ds[j] += 1
        memo2[cnt] = calc_cap(ds)
        print i, j, memo2[cnt], ds
        cnt+=1


    

    




1 1 0.000487618361305 [ 0.  2.  0.  0.  0.  0.  0.  0.  0.  0.]
1 2 0.0006919110623 [ 0.  1.  1.  0.  0.  0.  0.  0.  0.  0.]
1





    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-19-d92e3aab7241> in <module>()
      6         ds[i] += 1
      7         ds[j] += 1
----> 8         memo2[cnt] = calc_cap(ds)
      9         print i, j, memo2[cnt], ds
     10         cnt+=1

<ipython-input-7-73b756d7ce81> in calc_cap(ds)
      7         y[i] = tmp
      8     lin = linear_model.LinearRegression()
----> 9     lin.fit(b[T0:T0+T1,:], y[:T1])
     10     z = lin.predict(b[T0:,:])
     11 

/home/mattya/anaconda/lib/python2.7/site-packages/sklearn/linear_model/base.pyc in fit(self, X, y, n_jobs)
    297         else:
    298             self.coef_, self.residues_, self.rank_, self.singular_ = \
--> 299                 linalg.lstsq(X, y)
    300             self.coef_ = self.coef_.T
    301 

/home/mattya/anaconda/lib/python2.7/site-packages/scipy/linalg/basic.pyc in lstsq(a, b, cond, overwrite_a, overwrite_b)
    483         v, x, s, rank, work, info = gelss(
    484             a1, b1, cond=cond, lwork=lwork, overwrite_a=overwrite_a,
--> 485             overwrite_b=overwrite_b)
    486 
    487     else:

KeyboardInterrupt: 





    




 3 0.000670852447598 [ 0.  1.  0.  1.  0.  0.  0.  0.  0.  0.]

In [160]:

    

memo3 = zeros(1000, dtype=float)
cnt = 0
for i in range(0, 10):
    for j in range(i, 10):
        for k in range(j, 10):
            ds = zeros(10)
            ds[i] += 1
            ds[j] += 1
            ds[k] += 1
            
            memo3[cnt] = calc_cap(ds)
            print memo3[cnt], ds
            cnt+=1


    

    




0.00155398423335 [ 3.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
0.00134786499496 [ 2.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
0.00225959965066 [ 2.  0.  1.  0.  0.  0.  0.  0.  0.  0.]
0.00186634351455 [ 2.  0.  0.  1.  0.  0.  0.  0.  0.  0.]
0.00148023800149 [ 2.  0.  0.  0.  1.  0.  0.  0.  0.  0.]
0.00151986042376 [ 2.  0.  0.  0.  0.  1.  0.  0.  0.  0.]
0.00185838639106 [ 2.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
0.0013644368556 [ 2.  0.  0.  0.  0.  0.  0.  1.  0.  0.]
0.00189895934131 [ 2.  0.  0.  0.  0.  0.  0.  0.  1.  0.]
0.00227176924649 [ 2.  0.  0.  0.  0.  0.  0.  0.  0.  1.]
0.00155831606001 [ 1.  2.  0.  0.  0.  0.  0.  0.  0.  0.]
0.00149572447962 [ 1.  1.  1.  0.  0.  0.  0.  0.  0.  0.]
0.00166261458535 [ 1.  1.  0.  1.  0.  0.  0.  0.  0.  0.]
0.00110537581697 [ 1.  1.  0.  0.  1.  0.  0.  0.  0.  0.]
0.00160087484798 [ 1.  1.  0.  0.  0.  1.  0.  0.  0.  0.]
0.00179791788575 [ 1.  1.  0.  0.  0.  0.  1.  0.  0.  0.]
0.00127392963775 [ 1.  1.  0.  0.  0.  0.  0.  1.  0.  0.]
0.00244312364674 [ 1.  1.  0.  0.  0.  0.  0.  0.  1.  0.]
0.00135170833316 [ 1.  1.  0.  0.  0.  0.  0.  0.  0.  1.]
0.00130726729265 [ 1.  0.  2.  0.  0.  0.  0.  0.  0.  0.]
0.00135597096189 [ 1.  0.  1.  1.  0.  0.  0.  0.  0.  0.]
0.00190647758583 [ 1.  0.  1.  0.  1.  0.  0.  0.  0.  0.]
0.00221204280965 [ 1.  0.  1.  0.  0.  1.  0.  0.  0.  0.]
0.00175635553196





    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-160-3490d928b556> in <module>()
      9             ds[k] += 1
     10 
---> 11             memo3[cnt] = calc_cap(ds)
     12             print memo3[cnt], ds
     13             cnt+=1

<ipython-input-157-c3ac123852f0> in calc_cap(ds)
      3     for i in range(T1):
      4         tmp = 1.0
----> 5         for j in range(len(ds)):
      6             tmp *= legendre(ds[j], a[T0+i-j])
      7         y[i] = tmp

KeyboardInterrupt: 





    




 [ 1.  0.  1.  0.  0.  0.  1.  0.  0.  0.]

In [ ]:

    

print sum(memo[memo>0])
#print sum(memo2[memo2>0])
print sum(memo3[memo3>0])


    

In [42]:

    




    

    
Out[42]:





0.00037896847984852933

In [111]:

    

figure(10)


    

    
Out[111]:

In [112]:

    

yscale("log")


    

In [115]:

    

ylim(0.05, 1)


    

    
Out[115]:





(0.05, 1)

In [116]:

    

xscale("log")


    

In [117]:

    

xlim(10, 100)


    

    
Out[117]:





(10, 100)

In [ ]: