In [7]:
from numpy import *
from scipy import *
from sklearn import *

In [10]:
T0 = 100000
T1 = 100000

In [11]:
N = 50
alpha = 0.1
maxSVD = 0.9

In [51]:
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 [79]:
a = loadtxt("input.txt")
b = loadtxt("output.txt")

In [80]:
figure()
plot(b[:1000,0])
plot(b[:1000,1])
plot(b[:1000,2])
plot(b[:1000,3])


Out[80]:
[<matplotlib.lines.Line2D at 0x22ed44ec>]

In [40]:
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(b[T0:T0+T1,:], y[:T1])
    z = lin.predict(b[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 [81]:
memo = zeros(25, dtype=float)
for i in range(25):
    ds = zeros(i+1)
    ds[i] = 1
    memo[i] = calc_cap(ds)
    print i, memo[i]


0 0.999999999957
1 0.999999998892
2 0.999999988983
3 0.999999928036
4 0.999999251472
5 0.99999658378
6 0.999989120219
7 0.999931878074
8 0.999668880919
9 0.997681248958
10 0.981620387916
11 0.915700783508
12 0.793992485019
13 0.599603333437
14 0.236257532167
15 0.0533869184322
16 0.00771134606204
17 0.00308611539721
18 0.000923740776231
19 0.000644171881906
20 0.000430976736178
21 0.000425595121443
22 0.000423302495533
23 0.000575855906903
24 0.000510210990413

In [24]:
clf()
plot(memo)


Out[24]:
[<matplotlib.lines.Line2D at 0x1581720c>]

In [72]:
memo2 = zeros(100, dtype=float)
cnt = 0
for i in range(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


0 0 0.964576701958 [ 2.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
0 1 0.956209029382 [ 1.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
0 2 0.889166580813 [ 1.  0.  1.  0.  0.  0.  0.  0.  0.  0.]
0 3 0.636553683949 [ 1.  0.  0.  1.  0.  0.  0.  0.  0.  0.]
0 4 0.168601072056 [ 1.  0.  0.  0.  1.  0.  0.  0.  0.  0.]
0 5 0.0106816259871 [ 1.  0.  0.  0.  0.  1.  0.  0.  0.  0.]
0 6 0.00145794990133 [ 1.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
0 7 0.000604621221254 [ 1.  0.  0.  0.  0.  0.  0.  1.  0.  0.]
0 8 0.000434421254021 [ 1.  0.  0.  0.  0.  0.  0.  0.  1.  0.]
0 9 0.000626507604999 [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  1.]
1 1 0.959019250798 [ 0.  2.  0.  0.  0.  0.  0.  0.  0.  0.]
1 2 0.921782927423 [ 0.  1.  1.  0.  0.  0.  0.  0.  0.  0.]
1 3 0.663972546643 [ 0.  1.  0.  1.  0.  0.  0.  0.  0.  0.]
1 4 0.119709354627 [ 0.  1.  0.  0.  1.  0.  0.  0.  0.  0.]
1 5 0.0120620908073 [ 0.  1.  0.  0.  0.  1.  0.  0.  0.  0.]
1 6 0.00123649620396 [ 0.  1.  0.  0.  0.  0.  1.  0.  0.  0.]
1 7 0.000519191776608 [ 0.  1.  0.  0.  0.  0.  0.  1.  0.  0.]
1 8 0.000712736821787 [ 0.  1.  0.  0.  0.  0.  0.  0.  1.  0.]
1 9 0.000499730137409 [ 0.  1.  0.  0.  0.  0.  0.  0.  0.  1.]
2 2 0.889871498505 [ 0.  0.  2.  0.  0.  0.  0.  0.  0.  0.]
2 3 0.56882780048 [ 0.  0.  1.  1.  0.  0.  0.  0.  0.  0.]
2
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-72-7f44e09aac29> 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-40-929a6850de8a> in calc_cap(ds)
      4         tmp = 1.0
      5         for j in range(len(ds)):
----> 6             tmp *= legendre(ds[j], a[T0+i-j])
      7         y[i] = tmp
      8     lin = linear_model.LinearRegression()

<ipython-input-51-8e7ffb0a34f1> in legendre(n, x)
      1 def legendre(n, x):
----> 2     if n==0:
      3         return 1
      4     elif n==1:
      5         return x * (2*n+1)/2

KeyboardInterrupt: 
 4 0.0853456042745 [ 0.  0.  1.  0.  1.  0.  0.  0.  0.  0.]

In [68]:
memo3 = zeros(100, dtype=float)
cnt = 0
for i in range(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.998217339497 [ 3.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
0.999748887038 [ 2.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
0.996188485693 [ 2.  0.  1.  0.  0.  0.  0.  0.  0.  0.]
0.951696516942 [ 2.  0.  0.  1.  0.  0.  0.  0.  0.  0.]
0.827240250261 [ 2.  0.  0.  0.  1.  0.  0.  0.  0.  0.]
0.648086635414 [ 2.  0.  0.  0.  0.  1.  0.  0.  0.  0.]
0.288042627064 [ 2.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
0.257502185456 [ 2.  0.  0.  0.  0.  0.  0.  1.  0.  0.]
0.14916823708 [ 2.  0.  0.  0.  0.  0.  0.  0.  1.  0.]
0.0204884710279 [ 2.  0.  0.  0.  0.  0.  0.  0.  0.  1.]
0.995314651278 [ 1.  2.  0.  0.  0.  0.  0.  0.  0.  0.]
0.995880480688 [ 1.  1.  1.  0.  0.  0.  0.  0.  0.  0.]
0.919815115745 [ 1.  1.  0.  1.  0.  0.  0.  0.  0.  0.]
0.759317289821 [ 1.  1.  0.  0.  1.  0.  0.  0.  0.  0.]
0.648640321564 [ 1.  1.  0.  0.  0.  1.  0.  0.  0.  0.]
0.423205146712 [ 1.  1.  0.  0.  0.  0.  1.  0.  0.  0.]
0.10348933009 [ 1.  1.  0.  0.  0.  0.  0.  1.  0.  0.]
0.0524375029976 [ 1.  1.  0.  0.  0.  0.  0.  0.  1.  0.]
0.0172942935679 [ 1.  1.  0.  0.  0.  0.  0.  0.  0.  1.]
0.956911997112 [ 1.  0.  2.  0.  0.  0.  0.  0.  0.  0.]
0.832920371232 [ 1.  0.  1.  1.  0.  0.  0.  0.  0.  0.]
0.712155615974 [ 1.  0.  1.  0.  1.  0.  0.  0.  0.  0.]
0.257924888349 [ 1.  0.  1.  0.  0.  1.  0.  0.  0.  0.]
0.179586879116 [ 1.  0.  1.  0.  0.  0.  1.  0.  0.  0.]
0.0833628327284 [ 1.  0.  1.  0.  0.  0.  0.  1.  0.  0.]
0.0459238581594 [ 1.  0.  1.  0.  0.  0.  0.  0.  1.  0.]
0.012549531655 [ 1.  0.  1.  0.  0.  0.  0.  0.  0.  1.]
0.517941007694 [ 1.  0.  0.  2.  0.  0.  0.  0.  0.  0.]
0.419349097963 [ 1.  0.  0.  1.  1.  0.  0.  0.  0.  0.]
0.295563601135 [ 1.  0.  0.  1.  0.  1.  0.  0.  0.  0.]
0.0728154642586 [ 1.  0.  0.  1.  0.  0.  1.  0.  0.  0.]
0.0308358824707 [ 1.  0.  0.  1.  0.  0.  0.  1.  0.  0.]
0.0074974061144 [ 1.  0.  0.  1.  0.  0.  0.  0.  1.  0.]
0.00278273222253 [ 1.  0.  0.  1.  0.  0.  0.  0.  0.  1.]
0.055534897651 [ 1.  0.  0.  0.  2.  0.  0.  0.  0.  0.]
0.041320316726 [ 1.  0.  0.  0.  1.  1.  0.  0.  0.  0.]
0.0125770965192 [ 1.  0.  0.  0.  1.  0.  1.  0.  0.  0.]
0.00476450380907 [ 1.  0.  0.  0.  1.  0.  0.  1.  0.  0.]
0.000977614958047 [ 1.  0.  0.  0.  1.  0.  0.  0.  1.  0.]
0.000662700696793 [ 1.  0.  0.  0.  1.  0.  0.  0.  0.  1.]
0.00503552456096 [ 1.  0.  0.  0.  0.  2.  0.  0.  0.  0.]
0.00544973684295 [ 1.  0.  0.  0.  0.  1.  1.  0.  0.  0.]
0.00304721856355 [ 1.  0.  0.  0.  0.  1.  0.  1.  0.  0.]
0.000975081235223 [ 1.  0.  0.  0.  0.  1.  0.  0.  1.  0.]
0.000475841673749 [ 1.  0.  0.  0.  0.  1.  0.  0.  0.  1.]
0.0017933783349 [ 1.  0.  0.  0.  0.  0.  2.  0.  0.  0.]
0.00158334626013 [ 1.  0.  0.  0.  0.  0.  1.  1.  0.  0.]
0.00096260485173 [ 1.  0.  0.  0.  0.  0.  1.  0.  1.  0.]
0.000316752062274 [ 1.  0.  0.  0.  0.  0.  1.  0.  0.  1.]
0.000410561273596 [ 1.  0.  0.  0.  0.  0.  0.  2.  0.  0.]
0.000595627199714 [ 1.  0.  0.  0.  0.  0.  0.  1.  1.  0.]
0.000333049606517 [ 1.  0.  0.  0.  0.  0.  0.  1.  0.  1.]
0.000479649625765 [ 1.  0.  0.  0.  0.  0.  0.  0.  2.  0.]
0.000539081960675 [ 1.  0.  0.  0.  0.  0.  0.  0.  1.  1.]
0.000713934565276 [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  2.]
0.987476452906 [ 0.  3.  0.  0.  0.  0.  0.  0.  0.  0.]
0.996472464787 [ 0.  2.  1.  0.  0.  0.  0.  0.  0.  0.]
0.982886862356 [ 0.  2.  0.  1.  0.  0.  0.  0.  0.  0.]
0.765666924995 [ 0.  2.  0.  0.  1.  0.  0.  0.  0.  0.]
0.523273179158 [ 0.  2.  0.  0.  0.  1.  0.  0.  0.  0.]
0.321250179965
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-68-53538558cf34> in <module>()
      9             ds[k] += 1
     10 
---> 11             memo3[cnt] = calc_cap(ds)
     12             print memo3[cnt], ds
     13             cnt+=1

<ipython-input-40-929a6850de8a> in calc_cap(ds)
      4         tmp = 1.0
      5         for j in range(len(ds)):
----> 6             tmp *= legendre(ds[j], a[T0+i-j])
      7         y[i] = tmp
      8     lin = linear_model.LinearRegression()

KeyboardInterrupt: 
 [ 0.  2.  0.  0.  0.  0.  1.  0.  0.  0.]

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


19.5322570684
23.6183428599

In [42]:



Out[42]:
0.00037896847984852933

In [ ]: