memoring

import numpy as np
import scipy.signal as signal
import sys
from tools import unshape

NA = np.newaxis

def runave(rdata, length):
    """runnnig average
    
    Arguments:

       'array'  -- first dimension must be time. 

       'length' -- runnning mean length

    """
    ntim = rdata.shape[0]
    if ntim < length:
        print "input array first dimension length must be larger than length."
        sys.exit()        

    if rdata.ndim == 1:
        if length%2 != 0:
            w = np.ones(length)/float(length)
        else:
            w = np.r_[0.5, np.ones(length-1), 0.5]/float(length)                
        runarray = np.ma.array(signal.convolve(rdata, w, 'same'))
        runarray[:length/2] = np.ma.masked
        runarray[-length/2:] = np.ma.masked
    elif rdata.ndim == 2:
        if length%2 != 0:
            w = np.vstack((np.zeros(length),np.ones(length),np.zeros(length)))/float(length)
        else:
            w = np.vstack((np.zeros(length+1),np.r_[0.5, np.ones(length-1), 0.5],\
                           np.zeros(length+1)))/float(length)
        runarray = np.ma.array(signal.convolve2d(rdata, w, 'same'))
        runarray[:length/2,:] = np.ma.masked
        runarray[-length/2:,:] = np.ma.masked
    else :
        rdata, oldshape = unshape(rdata) 
        if length%2 != 0:
            w = np.vstack((np.zeros(length),np.ones(length),np.zeros(length)))/float(length)
        else:
            w = np.vstack((np.zeros(length+1),np.r_[0.5, np.ones(length-1), 0.5],\
                           np.zeros(length+1)))/float(length)
        runarray = np.ma.array(signal.convolve2d(rdata, w, 'same'))
        runarray[:length/2,:] = np.ma.masked
        runarray[-length/2:,:] = np.ma.masked
        runarray = runarray.reshape(oldshape)

    return runarray

ベクトル同士の線形回帰ならnumpy.linalg.lstsq や scipy.stats.linregress がよい。同じ時系列をもつ多次元配列をまとめて線形回帰するにはnumpy.polyfitが使える。

def regression(x, y):
    """
    regression y = ax + b
    """
    if x.ndim != 1:
        print "error x must be 1-D array"
        sys.exit()
    if y.ndim >= 2:
        y, oldshape = unshape(y)
        p = np.polyfit(x,y,1)
        a = p[0].reshape(oldshape[1:])
        b = p[1].reshape(oldshape[1:])
    else:
        p = np.polyfit(x,y,1)
        a = p[0]
        b = p[1]

    return a, b

データ行列を特異値分解することにより求める。数学的な原理はこちらを参考。こんな感じになるか（要確認）。

import numpy
import scipy.linalg
import math

NA = numpy.newaxis
PI = 4.*math.atan(1.)

#経度144、緯度37、時間31のGrADS形式のバイナリを例に
f = open("ファイルパス","rb")
data = numpy.fromfile(f, "<f").reshape(31,37,144)
f.close()

#31年平均を気候値として偏差を求める
residual = data - data.mean(axis=0)

#緯度に応じた面積重みをつける
lat = numpy.arange(0., 91., 2.5) 
input = residual*numpy.sqrt(numpy.cos(2.0*PI/360.0*lat))[NA,:,NA]
input = input.reshape(31,37*144)

#特異値分解
A, Lh, E = scipy.linalg.svd(input, full_matrices=False)

#固有値が並んだ配列
lambdas = Lh*Lh/(len(residual)-1)
#固有ベクトル[w1, w2, ... ]
EOFs = numpy.transpose(E)
#PC
PCs = A*Lh

#再規格化する場合
EOFs = EOFs * numpy.sqrt(lambdas)[NA,:]
PCs = PCs / numpy.sqrt(lambdas)[NA,:]

#緯度重みを考慮
EOFs = EOFs.reshape(37,144,len(lambdas)) / numpy.sqrt(numpy.cos(2.0*PI/360.0*lat))[NA,:,NA]

#寄与率
explained = lambdas / numpy.add.reduce(lambdas)

#第mモードのEOF、PC、寄与率はEOFs[:,:,m], PCs[:,m]、explained[m]になる（はず）

分散共分散行列の固有値問題を用いる場合は以下のような感じになるか。SVDを用いる方が早い気がする。