センター試験の群数列の問題をPythonで解いてみた。(2)

　前回の続きです。途中の式を最小二乗法を使って求めてみました。(^_^;

$\frac{1}{2}|\frac{1}{3},\frac{2}{3}|\frac{1}{4},\frac{2}{4},\frac{3}{4}|\frac{1}{5},\dots$
(2-1)数列{ $a_n$ }において $(\frac{1}{k})$ が初めて現れる項、第 $M_k$ 項
(2-2)数列{ $a_n$ }において $(\frac{k-1}{k})$ が初めて現れる項、第 $N_k$ 項
(3-1)数列{ $a_n$ }の第 $M_k$ 項から第 $N_k$ 項までの和 $S_k$
(3-2)数列{ $a_n$ }の初項から第 $N_k$ 項までの和 $T_k$

ただし、(3-1)及び(3-2)の和にそれぞれ名前を付けました。
　たとえば、第 $M_k$ 項の場合、式の形が与えられているので、
　 $M(k)=pk^2+qk+r$
とおいて、実際のデータから連立方程式を作って解いてしまえば難しいことは考えなくても機械的に解けてしまいます。(^_^;

n  1     2     3     4     5     6     7        
   1  |  1     2  |  1     2     3  |  1        
  --- | --- , --- | --- , --- , --- | --- , ... 
   2  |  3     3  |  4     4     4  |  5        
 M(2)  M(3)        M(4)              M(5)       
 N(2)        N(3)              N(4)

k	2	3	4
Mk	1	2	4
Nk	1	3	6
Sk	0.5	1.0	1.5
Tk	0.5	1.5	3.0

　よって、データから、
　 $M(2)= 4p+2q+r=1$
　 $M(3)= 9p+3q+r=2$
　 $M(4)=16p+4q+r=4$
　これを解けば答えは出ますが、式を作るのが面倒なので、データだけ与えれば解けてしまう最小二乗法を使うことにしました。(^_^;
　ここで、フィットさせたい関数は2次式と1次式ですが、まとめて2次式にしておけば、2次の係数が0のときは1次式にもなります。
　ちなみに、正攻法で解くとすれば、第ｉ群の項数は、ｉ[個]で、群数列の第ｉ群のｊ番目の項が原数列の第ｎ項だとすると、次式が成り立ちます。
　 $n=\frac{i(i-1)}{2}+j$
　これとk=i+1∴i=k-1から、
　j=1　のとき、 $M_k=\frac{(k-1)(k-2)}{2}+1=\frac{k(k-3)}{2}+2$
　j=k-1のとき、 $N_k=\frac{(k-1)(k-2)}{2}+k-1=\frac{k(k-1)}{2}$
　また、数列{ $a_n$ }の第 $M_k$ 項から第 $N_k$ 項までの和は、
　 $\frac{1}{k}+\frac{2}{k}+\dots+\frac{k-1}{k}=\sum_{j=1}^{k-1}\frac{j}{k}=\frac{1}{k}\sum_{j=1}^{k-1}j=\frac{1}{2}(k-1)$
　よって、数列{ $a_n$ }の初項から第 $N_k$ 項までの和は、
　 $\sum_{i=2}^{k}\frac{1}{2}(i-1)=\frac{1}{2}\left(1+2+\dots+k-1\right)=\frac{1}{2}\sum_{i=1}^{k-1}i=\frac{1}{4}k(k-1)$

P.S.
　初め、連立方程式で解くと式を立てるのが面倒かと思っていましたが、同じ2次式で統一すれば、まとめて扱えるようなので検算用としてプログラムに追加しておきました。(^_^;

● GrpSeqCt2.py

# coding: UTF-8
# GrpSeqCt2.py

import numpy as np
from scipy.optimize import leastsq
from fractions import Fraction

def fit_fnc(param,x):
    a,b,c = param
    return a*x*x+b*x+c

# residual
def rsd_fnc(param,x,y):
    return y-fit_fnc(param,x)

def main():
    I = 3   # 3以上の整数
    a = ['0/1']
    M,N,S,T = [],[],[],[]
    n = t = 0
    for i in range(1,I+1):
        s = 0
        for j in range(1,i+1):
            frc = str(j)+'/'+str(i+1)
            a.append(frc)
            s+= Fraction(frc)
            n+= 1
            if j==1: M.append(n)
            if j==i: N.append(n)
        S.append(float(s))
        t+=s
        T.append(float(t))
##    print("%s\n%s\n%s\n%s\n%s\n"%(a,M,N,S,T))
    Pk = np.array(range(2,I+2))
    Pm,Pn = np.array(M),np.array(N)
    Ps,Pt = np.array(S),np.array(T)
    param0 = [0.,0.,0.]     # a,b,cの初期値
    Mk = leastsq(rsd_fnc,param0,args=(Pk,Pm))[0]
    Nk = leastsq(rsd_fnc,param0,args=(Pk,Pn))[0]
    Sk = leastsq(rsd_fnc,param0,args=(Pk,Ps))[0]
    Tk = leastsq(rsd_fnc,param0,args=(Pk,Pt))[0]
    print('(2-1) Mk=(% f)k^2+(% f)k+(% f)'%(Mk[0],Mk[1],Mk[2]))
    print('(2-2) Nk=(% f)k^2+(% f)k+(% f)'%(Nk[0],Nk[1],Nk[2]))
    print('(3-1) Sk=(% f)k^2+(% f)k+(% f)'%(Sk[0],Sk[1],Sk[2]))
    print('(3-2) Tk=(% f)k^2+(% f)k+(% f)'%(Tk[0],Tk[1],Tk[2]))

    print(u'//--- 検算 ---//')
    A = np.array([[ 4, 2, 1],
                  [ 9, 3, 1],
                  [16, 4, 1]])
    bM,bN,bS,bT = Pm[:3],Pn[:3],Ps[:3],Pt[:3]
    Mk = np.linalg.solve(A,bM)
    Nk = np.linalg.solve(A,bN)
    Sk = np.linalg.solve(A,bS)
    Tk = np.linalg.solve(A,bT)
    print('(2-1) Mk=(% f)k^2+(% f)k+(% f)'%(Mk[0],Mk[1],Mk[2]))
    print('(2-2) Nk=(% f)k^2+(% f)k+(% f)'%(Nk[0],Nk[1],Nk[2]))
    print('(3-1) Sk=(% f)k^2+(% f)k+(% f)'%(Sk[0],Sk[1],Sk[2]))
    print('(3-2) Tk=(% f)k^2+(% f)k+(% f)'%(Tk[0],Tk[1],Tk[2]))

if __name__ == '__main__':
    main()

●実行結果

(2-1) Mk=( 0.500000)k^2+(-1.500000)k+( 2.000000)
(2-2) Nk=( 0.500000)k^2+(-0.500000)k+(-0.000000)
(3-1) Sk=(-0.000000)k^2+( 0.500000)k+(-0.500000)
(3-2) Tk=( 0.250000)k^2+(-0.250000)k+( 0.000000)
//--- 検算 ---//
(2-1) Mk=( 0.500000)k^2+(-1.500000)k+( 2.000000)
(2-2) Nk=( 0.500000)k^2+(-0.500000)k+(-0.000000)
(3-1) Sk=( 0.000000)k^2+( 0.500000)k+(-0.500000)
(3-2) Tk=( 0.250000)k^2+(-0.250000)k+(-0.000000)

※参考URL
●大学入試センター試験｜解答速報2016｜予備校の東進
●センター試験の群数列の問題をPythonで解いてみた。
●最小二乗法 - Python Project - Seesaa Wiki（ウィキ）
●pythonでフィッティングをする - おっぱいそん！
●scipyで最小二乗法 - daharuの日記
●Python Numpy Scipy サンプルコード: 線形連立方程式, LU 分解 | org-技術