有哪些向量化寫法讓你拍案叫絕？
在很多數據分析中（還有函數式編程中），向量式（vectorized）寫法可以大大提升計算效率（比如R、Python Pandas）裡面，而且向量式寫法也更加容易並行運算。有哪些你寫過或者見過的向量式寫法令你拍案叫絕的？我自己印象深刻的是向量式統計連續遞增遞減序列
我來裝回b吧。以下代碼出自這 PRML/PRMLT.在老版本Matlab上，基本都比內置函數快個1~2個數量級，新版Matlab沒測過了。兩組數據歐式距離（閻小鵬上面提了）
function D=sqdist(X, Y)
D=bsxfun(@plus,dot(X,X,1)",dot(Y,Y,1))-2*(X"*Y);
function [label, energy, m] = litekmeans(X, k)
% Perform k-means clustering.
% Input:
% X: d x n data matrix
% init: k number of clusters or label (1 x n vector)
% Output:
% label: 1 x n cluster label
% energy: optimization target value
% m: d x k center of clusters
% Written by Mo Chen (sth4nth@gmail.com).
n = size(X,2);
idx = 1:n;
last = zeros(1,n);
label = ceil(k*rand(1,n)); % random initialization
while any(label ~= last)
[~,~,last(:)] = unique(label); % remove empty clusters
E = sparse(idx,last,1); % transform label into indicator matrix
m = X*(E./sum(E,1)); % compute centers
[val,label] = min(dot(m,m,1)"/2-m"*X,[],1); % assign labels
end
energy = dot(X(:),X(:),1)+2*sum(val);
function [label, energy, model] = knKmeans(X, k, kn)
n = size(X,2);
K = kn(X,X);
last = zeros(1,n);
label = ceil(k*rand(1,n));
while any(label ~= last)
[~,~,last(:)] = unique(label); % remove empty clusters
E = sparse(last,1:n,1);
E = E./sum(E,2);
T = E*K;
[val, label] = max(T-dot(T,E,2)/2,[],1);
end
energy = trace(K)-2*sum(val);
function [label, energy，index] = kmedoids(X, init)
X = X-mean(X,2); % reduce chance of numerical problems
v = dot(X,X,1);
D = v+v"-2*(X"*X); % Euclidean distance matrix
D(sub2ind([d,d],1:d,1:d)) = 0; % reduce chance of numerical problems
last = zeros(1,n);
label = ceil(k*rand(1,n));
while any(label ~= last)
[~,~,last(:)] = unique(label); % remove empty clusters
[~, index] = min(D*sparse(1:n,last,1),[],1); % find k medoids
[val, label] = min(D(index,:),[],1); % assign labels
end
energy = sum(val);
function [Q, R] = gson(X)
[d,n] = size(X);
m = min(d,n);
R = zeros(m,n);
Q = zeros(d,m);
for i = 1:m
R(1:i-1,i) = Q(:,1:i-1)"*X(:,i);
v = X(:,i)-Q(:,1:i-1)*R(1:i-1,i);
R(i,i) = norm(v);
Q(:,i) = v/R(i,i);
end
R(:,m+1:n) = Q"*X(:,m+1:n);
function y = loggausspdf(X, mu, Sigma)
d = size(X,1);
X = bsxfun(@minus,X,mu);
[U,p]= chol(Sigma);
if p ~= 0
error("ERROR: Sigma is not PD.");
end
Q = U"X;
q = dot(Q,Q,1); % quadratic term (M distance)
c = d*log(2*pi)+2*sum(log(diag(U))); % normalization constant
y = -(c+q)/2;
function [u, w] = nbBern(X, t)
n = size(X,2);
E = sparse(1:n,t,1);
nk = sum(E,1);
w = full(nk/n);
mu = X*(E./nk);
function [z, llh] = hmmViterbi(x, A, E, s)
k = numel(s);
Z = zeros(k,n);
Z(:,1) = 1:k;
v = s(:)+M(:,1);
for t = 2:n
[v,idx] = max(bsxfun(@plus,A,v),[],1); % 13.68
v = v(:)+M(:,t);
Z = Z(idx,:);
Z(:,t) = 1:k;
end
[llh,idx] = max(v);
z = Z(idx,:);
function [gamma, alpha, beta, c] = hmmSmoother(M, A, s)
[K,T] = size(M);
At = A";
c = zeros(1,T);
alpha = zeros(K,T);
[alpha(:,1),c(1)] = normalize(s.*M(:,1),1);
for t = 2:T
[alpha(:,t),c(t)] = normalize((At*alpha(:,t-1)).*M(:,t),1); % 13.59
end
beta = ones(K,T);
for t = T-1:-1:1
beta(:,t) = A*(beta(:,t+1).*M(:,t+1))/c(t+1); % 13.62
end
gamma = alpha.*beta; % 13.64
function label = litekmeans(X, k)
% Perform k-means clustering.
% X: d x n data matrix
% k: number of seeds
% Written by Michael Chen (sth4nth@gmail.com).
n = size(X,2);
last = 0;
label = ceil(k*rand(1,n)); % random initialization
while any(label ~= last)
[u,~,label] = unique(label); % remove empty clusters
k = length(u);
E = sparse(1:n,label,1,n,k,n); % transform label into indicator matrix
m = X*(E*spdiags(1./sum(E,1)",0,k,k)); % compute m of each cluster
last = label;
[~,label] = max(bsxfun(@minus,m"*X,dot(m,m,1)"/2),[],1); % assign samples to the nearest centers
end
[~,~,label] = unique(label);
這個Kmeans的Matlab代碼讓我印象深刻，比matlab內置的kmeans函數快一個數量級。
最新版本在這裡Kmeans Clustering關於全排列（Permutation）問題，@親愛的龍哥 龍哥給出的演算法確實強大。沒去看的知友可以去先了解下。
有哪些向量化寫法讓你拍案叫絕？ - 親愛的龍哥的回答 - 知乎我用遞歸稍微改寫了龍哥的演算法，主要是為了方便大家理解龍哥演算法的核心思想。#R code
#這是龍哥的演算法
permns_longge &<- function(n,vec = 1:n)
{
ans = matrix(1, nrow = factorial(n), ncol = n)
faci = 1
if(n&>1)
{
for(i in 2:n)
{
faci_1 = faci
faci = faci_1 * i
series = 1:i
head = faci - faci_1 + 1
for(p in i:1)
{
rvec = series[-p]
ans[head:(head + faci_1 - 1),2:i] = rvec[ans[1:faci_1,1:(i-1)]]
head = head - faci_1
}
ans[1:faci,1] = rep(1:i, each = faci_1)
}
}
ans = matrix(vec[ans], ncol = n)
ans
}
#這是我用遞歸修改後的
permns &<- function(n,vec = 1:n)
{
if(n==1)return(vec)
if(n&>1)
{
facn_1 = factorial(n-1) #factorial(n-1)
facn = facn_1 * n #factorial(n)
ans = matrix(1, nrow = facn, ncol = n)
series = 1:n
head = facn - facn_1 + 1
permnsn_1 = permns(n-1) #the (N-1)th result
for(p in n:1)
{
rvec = series[-p]
ans[head:(head + facn_1 - 1),2:n] = rvec[permnsn_1]
head = head - facn_1
}
ans[1:facn,1] = rep(1:n, each = facn_1)
ans = matrix(vec[ans], ncol = n)
ans
}
}
使用遞歸的話，因為沒有嵌套循環所以方便大家理解啦~實際上就是把1到n-1的全排列結果拿來填充1到n的全排列矩陣的後面n-1列，再把第一列的1到n補上
比較如下：&> system.time(permns_longge(10))
用戶 系統 流逝
1.33 0.30 1.63
&> system.time(permns(10))
用戶 系統 流逝
0.91 0.17 1.08
計算1到10的全排列，在我的PC（i7-4720HQ CPU 12G RAM）上運行兩種方法，都是1秒多，遞歸改寫的平均能快0.3到0.5秒吧。嗯，好像有那麼點用(*////▽////*)。
所以還是為了方便大家理解龍哥的精華。嗯，就是這樣。。。
另外，輸出全排列不只是對數字有效哦，例如：&> permns(3,c("a","b","c"))
[,1] [,2] [,3]
[1,] "a" "b" "c"
[2,] "a" "c" "b"
[3,] "b" "a" "c"
[4,] "b" "c" "a"
[5,] "c" "a" "b"
[6,] "c" "b" "a"
更新：補充了一下內容
對於1：上三角矩陣的形式；對於3： 關於Quant 通常需要什麼樣的教育背景和知識結構？ - 袁浩瀚的回答 統計第一題均勻更新過程次數期望隨著和的增加而收斂至本質上，矩陣化就是利用了流行科學計算語言里的批量生成功能，用內存換取時間。把複雜計算問題表示成矩陣運算問題，再利里科學計算語言包里已經優化過的基礎矩陣運算功能大幅減少計算時間。 這樣做的額外好處是，對數學工作者，該類語句都是具體清晰易懂的簡單表達式，執行的每一步都是一個矩陣運算，debug容易（數值上）。
一旦把所有問題轉化為矩陣語言，在有矩陣函數的前提下進入矩陣化編程領域，傳統的編程思維就可以統統扔一邊，全部轉化成計算題。對！在非數學工作程序員眼裡，我們都是邪教！我們就是要把這種邪教發揚光大！1.MC求期望或者敏感度時候的「二向箔」：對於一個基於一個利率模型的固收價格模擬我們總能找到三個維度：
補充 zero bond（貼現率）一次模擬的 三角形式：
上面這個矩陣可以很輕易地通過一些 諸如 tril 或者 rot90 的函數完成
這個時候，我們如果多建一個維度（也就是mc軸），再利用所有矩陣語言都可以修改對求和（SUM）均值（AVerage）等線性函數在矩陣里的方向（比如Python就可以用0或1或n，來決定對哪個維度使用該函數）。 我們求一些floating leg或者該過程下Z的期望價格的時候，完全可以利用對某個方向的「壓縮」，實現模擬的 「二向箔」：
2.做數據篩選處理的時候，巧用信號類函數的矩陣形式：
一般而言，矩陣語言里都有 選出某滿足某一條件元素 index的功能，但是具體通過這些條件改變矩陣的實現會有一些困難。但是妙用一些信號函數（階躍，衝激）會使這個過程簡單很多，比如如下代碼節選：
M&<-matrix(rnorm(MC,mean=0,sd=1))Mcube&<-kronecker(matrix(1,numlie,1),t(M))dim(Mcube)&<-c(numlie,MC)Mcube&<-kronecker(matrix(1,1,numcomp),Mcube)dim(Mcube)&<-c(numlie,numcomp,MC)rouMpart&<-roucube*Mcube
sqrtroupart&<-(1-roucube^2)^0.5compmian&<-matrix(rnorm(numcomp*MC,mean = 0,sd=1),nrow = MC,ncol = numcomp) ##a whole matrix of randomness,MC*numcomp
compcube&<-kronecker(compmian,matrix(1,numlie,1))
dim(compcube)&<-c(numlie,numcomp,MC)compcube[3,,] ##supposed to be all different random numbers##company value done
compvalue&<-rouMpart*Mcube+sqrtroupart*compcubedim(compvalue)##threshold cube
threshold&<-matrix(seq(0.01,0.1,0.01),nrow = numthreshold,ncol=1) ## [th1,..,th10] sequence, 10 in totalthresholdlie&<-kronecker(matrix(1,numrou,1),threshold) ## 10 [th1,..,th10] sequence, 100 in totalthresholdcube&<-kronecker(matrix(1,MC,numcomp),thresholdlie)dim(thresholdcube)&<-c(numlie,numcomp,MC)
thresholdcube[21:30,,1]##result cube
differencecube&<-compvalue-qnorm(thresholdcube)##count
H&<-Heaviside(-differencecube,a=0) ##less than zero will be 1, bigger will be zerocountDefaults&<-function(companyValues) {
# the function count default payoff calculates how many defaults you have # per one path of the 100 companies. #step 1 - filter the elements that are below the threshold #step 2 - count the length of the vector result&<-companyValues[companyValues==1] return(length(result)) #counts the number of elements below the threshold i.e. the number of defaults}mcresult&<-apply(H,c(1,3),countDefaults)mcresult ##add the company list into the roulie*MC matrixdim(mcresult)&<-c(numlie,MC)這段R代碼節選是實現copula下信用違約模擬的， 用到了1裡面二相箔，而且立方體的構造使用了張量積；然而精髓是黑體：用 階躍信號函數（heaviside）來刪選 違約（quantile大於某一個閾值的模擬值）， 這個閾值可以直接通過調整階躍信號的參數來實現，上面的例子我選擇了做減法取0點然後套用標準階躍信號；使用躍遷函數直接把數值編程0或1之後，直接加總即可求得違約數 。由於我的躍遷函數可以只針對一個維度或者兩個維度設計，所以最後的違約數也可以不應複雜的find或者條件式也可以兩句話得到。
同時，如果想進行滿足條件篩選，也可以構造衝激信號（dirac ）來把不要的信息都過濾掉，留下來的又是一個0/1「計數」！。
3.使用累加和Unique函數簡化累積或障礙衍生品的運算：這個是我曾經試圖用數值法求解 @袁浩瀚 袁老師20題里 Quant 通常需要什麼樣的教育背景和知識結構？ - 袁浩瀚的回答統計第一題 時誤打誤撞想出來的我們構造了一個獨立增量過程random walk（隨便什麼比如泊松，二項）的時候，可以通過累加（matlab 裡面是 cumsum） 增量，模擬累積情況。當其觸線（比如設定一個障礙100）的時候，利用find和unique函數，先找到停時的第一個值，然後把這些滿足這些值得第一個調出來，即可簡便障礙或者累積的實現，這在非矩陣語言里，演算法實現非常繁瑣的
這個估計就是袁老師問題注釋里寫的：向量式統計連續遞增遞減序列， 吧
具體演算法示意圖：
補充：更新均勻更新過程隨著和增大，次數期望收斂的狀況,可以看到均勻的收斂的非常快，10以內就快0了
除此之外還有很多演算法上的壯舉是用矩陣完成的，包括： FFT 的遞歸實現， 線性微分運算元（對矩陣求導）的實現， 優化問題里的條件劃歸（凸優化，線性規劃 裡面都可以見到矩陣劃歸的影子），各種線性代數分解的矩陣實現法，還有 http://numerics.mathdotnet.com/t里的矩陣 lambda表達式賦值。
矩陣化的強大和魔性讓喵我深深著迷，每次寫代碼，別人看見都會說「黑貓又在畫方塊啊」 我覺得這就是矩陣化的魅力所在，噁心臃腫的循環結構，煩人的各種演算法，變成一個個立體形象的n維方塊。
想像黑客帝國的原名叫什麼： Matrices！ 給他一個矩陣，他能毀滅世界，VECTORIZATION！1. Permutation
我的思路在於數學歸納法.記perm_yifan &<- function(n){
n = floor(n);
if (n&<=1) stop("error");
if (n==3) {
return(matrix(c(3,3,1,2,1,2,1,2,3,3,2,1,2,1,2,1,3,3),ncol=3));
}
if (n==2) {
return(matrix(c(1,2,2,1),ncol=2));
}
tmp &<- perm_yifan(n-1);
re&<- matrix(0,n*nrow(tmp),ncol(tmp)+1)
ttmmpp &<- matrix(0,nrow(tmp),ncol(tmp)+1);
for (i in 1:n){
ttmmpp[,-i] = tmp;
ttmmpp[,i] = n
re[(1+(i-1)*nrow(tmp)):(i*nrow(tmp)),] = ttmmpp
}
re
}
對比 @GuQuQ 和 @親愛的龍哥 的演算法
&> system.time(perm_yifan(10))
user system elapsed
1.33 0.14 1.46
&> system.time(permns(10))
user system elapsed
1.70 0.33 2.03
運行環境Intel 5930K 64GB, Windows 10 Microsoft R Open
2. repmat問題
問題還是見 @親愛的龍哥 我評論裡面提到了要用R的特性, 這個坑是這麼填的. 我沒有認真優化那個中間變數的部分,因為應該不需要. 注意,我的BLAS是MKL. 默認的BLAS會更慢.
repmat_vect &<- function(A,m,n){ # my version
mm &<- nrow(A);
nn &<- ncol(A);
rep(A,n) -&> B
dim(B) &<- c(mm,n*nn)
rep(t(B),m) -&> BB
dim(BB) &<- c(m*mm,n*nn)
t(BB)
}
repmat_vect2 &<- function(X,m,n){ # stackoverflow: 19590541
mx = dim(X)[1]
nx = dim(X)[2]
matrix(t(matrix(as.double(X),mx,nx*n)),mx*m,nx*n,byrow=T)
}
repmat_alg &<- function(A,m,n){
matrix(1,m,n)%x%A
}
UU &<- matrix(rnorm(1000000),1000,1000)
&> system.time(re1 &<- repmat_vect(UU,10,10))
user system elapsed
2.00 0.26 2.28
&> system.time(re2 &<- repmat_vect2(UU,10,10))
user system elapsed
1.94 0.11 2.07
&> system.time(re3 &<- repmat_alg(UU,10,10))
user system elapsed
2.46 1.29 1.64
這題可以答上好一陣子了。。。
1 Permutation問題問題描述：如果給你1,2,3,...,9，讓你寫一個函數返回他們的9!=362880個全排列你會怎麼做？我當時在使用R語言開發一個包， 因為很多函數都需要生成全排列，而我又沒在網上找到一個合適的代碼，就只好自己寫了。不小心寫了一個R裡面最快的出來。我相信肯定有人會優化的比我好，但是目前開源的R代碼中我還沒有找到，可能這個問題太小眾了， 如果有，請告訴我，多謝。實現相同功能的函數， 我找到了如下幾個：combinat包的permn函數， gtools包的permutations函數，permute包的allPerms函數。 我的代碼如下：
我來解釋下這段代碼
這個問題最直觀的想法就是用遞歸， 你去Google搜索演算法，十有八九都是用遞歸的。在generate n個數的全排列的時候，我們會調用generate n-1個數的全排列，如此往複。遞歸的方法雖然直觀，但是當然是最慢的。如果我沒記錯的話，combinat包的permn函數和gtools的permutations函數都是用了遞歸的方法，不過剛才我查看了下源代碼似乎作者已經優化了。下面就有人給出來一個遞歸實現的例子（python）Algorithm to generate all possible permutations of a list?
如果我們把遞歸變成iterative的情況， 自然就優化了很多。
下面這個網站用了很多種語言來實現， 喜歡寫JAVA和C++的可以自行對比比遞歸方法能快多少。https://www.nayuki.io/page/next-lexicographical-permutation-algorithm
其實早在5年前的一次作業里， 我就用了一次向量化的方法來生成全排列， 想法如下：
假設我已經有了2個數的全排列矩陣：[1， 2; 2， 1]當我要生產3個數的全排列時，直接新建一個之全部為3的向量，並把他「夾在」兩個A之間，就變成了這樣：[ 1, 2, 3, 1, 2; 2, 1, 3, 2, 1]然後我對這個矩陣「切三刀」， 分別他的(1, 2, 3)列， （2， 3， 4）列， （3， 4， 5）列三個子矩陣：[1, 2, 3; 2, 1, 3][2, 3, 1;
1, 3, 2][3, 1, 2;
3 , 2, 1]OK了， 三個子矩陣縱向排在一起就是3元素的全排列了。 這個程序不難實現， 而且基本都是在矩陣的拼接和取子矩陣的操作。但是假如我們有更高的要求，我希望返回的全排列是 lexicographical的， 也就是說是按照從大到小的順序排列的， 那麼個這個演算法就需要做點修改了。（註：按照順序的意思是要按照（1,2,3) (1,3,2) (2,1,3) (2,3,1) (3,1,2) (3,2,1)這個順序來排列）於是乎我又有了新的想法。還是剛才的例子，假如我們已經生成了兩個元素的全排列， 我們要生成3個元素的全排列並且是按照順序來排列的。
要知道，這裡的1， 2並不僅僅是元素1， 2， 更重要的是他表示的是下標。 比如我要對兩個元素 v = c("tom", "mary")來實現全排列，那麼矩陣是應該是這個樣子的[v[1]， v[2]; v[2]， v[1]]也就是["tom","mary"; "mary", "tom"]有了這個認識， 那麼就OK了。 我們可以按照這個矩陣分別排列（2， 3） （1， 3） （1， 2）等二元組的全排列（注意我的代碼里的series[-p]就是在生成這個二元組）於是你的矩陣變成了這個樣子：[2， 3; 3， 2][1， 3; 3， 1][1， 2; 2， 1]我們只需要在他們的前面分別補充上缺失的第三個元素，1,1,2,2,3,3即可。 （代碼的ans[1:faci, 1] = rep(1:i, each = faci_1)就是在補充第一列）[1， 2， 3; 1， 3， 2][2， 1， 3; 2， 3， 1][3， 1， 2; 3， 2， 1]大功告成！2 repmat問題
這個演算法不是我原創，但是絕逼是能用上的最魔性的向量化演算法。眾所周知Matlab裡面有個repmat函數，相比R裡面的rep函數，repmat可以直接將矩陣複製M*N次函數如下：
Given a list of points, calculate the pairwise distances:
function dis = pdist3(X)
% input: X is a column vector of points
% return: dis is the pairwise distance matrix
X2 = sum(X.*X,2);
dis = sqrt(bsxfun(@plus, X2, X2") - 2*X*X");
end
比自帶的pdist不知道高到哪裡去了，速度和內存消耗都更優。
在計算分類正確率時：
mean((yhat ~= yte))
1. 在其他代數結構上的矩陣
正常的矩陣裡面的元素都是來自於域(現在的答案里幾乎只有實數域). 改為其他的代數結構比如某些特殊的semiring, 可以得到很多其他有用的結果.比如Gauss-Jordan-Floyd-Warshall-McNaughton-Yamada algorithm.
一個演算法可以同時解線性方程, 算圖上的最短距離, 求有限自動機的正則表達式.另一個類似的Fun with Semirings的文章paper: https://www.cl.cam.ac.uk/~sd601/papers/semirings.pdf 以及 slides: http://www.cl.cam.ac.uk/~sd601/papers/semirings-slides.pdf更多興趣有本書
Graphs, Dioids and Semirings2. 換個array processing的語言 (比如APL語言)
用APL語言會強迫自己一切都寫成向量的形式(以及一些古怪的APL的operator). K,J這些語言就發源於APL. 這可是讓Iverson拿到了圖靈獎的語言. 如今已經沒什麼人用了. 但是現在寫的話還是會培養向量形式的思維.
比如edit distance用APL可以這樣寫.
寫了一套回測平台，從指標計算到觸發到計算收益全是向量運算，速度確實快，不過比較耗內存
分組計算統計量
len = 1000; val = rand(len,1); group = randi(5,len,1);
group_summary = accumarray(group(:),val,[],@sum)
Floyd 演算法：
for k = 1:size(A,1)
A = min(A,A(:,k)+A(k,:));
end
1. 替換數據中的nan為平均值
some_raw = importdata("some.data", " ")";
[features, length] = size(some_raw);
navg = nanmean(some_raw);
nidx = isnan(some_raw);
midx = bsxfun(@and, nidx, sum(nidx) &< features);
some_raw(midx) = navg(nonzeros(bsxfun(@times, midx, 1:numel(navg))));
some_raw = some_raw";
B = bsxfun( @plus, ...
bsxfun( @minus, ...
bsxfun( @minus, D, sum(D, 1) / n ), ...
sum(D, 2) / n ), ...
sum(D(:)) / (n ^ 2) ) * (-0.5);
計算裡面常用函數之一的 diff() 如果用向量化表示可以為：
vector = seq(1,1000,1)
today = append(vector[1:999],NaN)
tommorow = append(NaN,vector[2:1000])
diff_result = tommorow - today
這個函數在做差分的時候經常用到，比如時間序列。
在實現最小二乘法做n次多項式擬合的時候，有這麼一段我認為寫得比較好的向量化代碼，實驗結果擬合速度甚至比matlab自帶的polyfit函數還快一些。
function coef=fitt(x, y, n)
m = length(x);
A = repmat(x", 1, 2*n);
c = repmat(1:2*n, m, 1);
vector =[m, sum(A.^c)];
index = meshgrid(1:n+1);
index = bsxfun(@plus,index, [0:n]");
B =[ones(1,m); repmat(x, n,1)];
cc = repmat([1,1:n]", 1, m);
b = sum(B.^cc * y", 2);
coef = fliplr((vector(index))");
上述函數完成最小二乘法演算法計算後為了保持與matlab自帶函數一樣的輸出格式，還另外將計算結果經裝置與filplr函數處理。
亮點在於在構造方程組的係數矩陣時並沒有直接用循環來構造每一個元素，而是生成了一個很有規律的係數向量，再向量化生成索引矩陣來最後生成係數矩陣。
與polyfit函數對比：x = linspace(0,1,5);
y = 1./(1+x);
tic
p1 = polyfit(x,y,4)
toc
tic
p2 = fitt(x,y,4)
toc
p1 =
0.1524 -0.5333 0.8667 -0.9857 1.0000
時間已過 0.003473 秒。
p2 =
0.1524 -0.5333 0.8667 -0.9857 1.0000
時間已過 0.000757 秒。
&>&>
pandas apply/ rolling_apply/expanding_apply/groupby……
我寫的策略基本都是用這幾個，再嵌入自己寫的function。
另外，第一次知道pandas.pct_change的時候確實驚艷了一下
不過pandas在對location(index/colume)的vectorize上沒有像對value那麼友善，還是至少要對行或列進行遍歷。也可能是我還沒找到正確的方法，有知道的大神求賜教啊！這裡簡單介紹下隱式馬爾可夫模型
閱讀原文：http://club.jr.jd.com/quant/topic/1080077
京東金融量化交流群：456448095
一．從文藝復興科技，大獎章和HMM說起
Renaissance Technologies（Ren Tech）這家對沖基金的名字在量化圈算是如雷貫耳了。創始人數學家James Simons帶領著一幫數學家、物理學家與市場博弈，其中最賺錢的投資組合就是1988年創立的Medallion Fund。1994年到2004年中期的年化收益高達71.8%，在全球金融危機的2008年，大部分對沖基金都虧損，而大獎章的return高達98.2%。
由這些天才般的數學家和物理學家以及一些超強的交易員構成的文藝復興科技和他們創建的神秘的大獎章基金，不禁讓外界猜測，他們究竟是搭建了一個怎樣的量化模型，才能用非絕對投機的方法戰勝市場呢？而隱式馬爾可夫模型（HMM）也由於一些原因被認為是他們最有可能運用的一個模型，也是本文即將介紹的最為精彩量化模型。
Ren Tech成立初期的創始人中有一位James的好朋友Lenny Baum，此人正是發明廣泛應用在語音識別等領域的Baum Welch演算法的那個Baum，演算法是用來確定隱式馬爾可夫模型中未知變數可能出現的概率，到今天可以說是廣泛應用於語言識別和金融領域了。並且，1993年加盟復興技術的劍橋大學數學博士尼可·帕特森就是全球HMM領域公認的專家。另外，James還雇了很多曾在IBM從事語音識別和自然語言處理的科學家來Ren Tech工作。以這些人的機器學習和文字信息處理的功底，很難不引人猜測Ren Tech很大可能就是利用了HMM或者HMM與其它統計方法如SVM結合來進行短期市場預測從而高頻交易的。《解密復興科技：基於隱蔽馬爾科夫模型的時序分析方法》 一書中介紹了這個Jim很有可能採用的演算法，含有詳細的公式推導和運算，並且附上了一些實證結果，可以參考看看。
二．認識模型
HMM與馬爾可夫模型相比，不同的地方在於隱藏變數。馬爾可夫鏈是具有有限記憶屬性的模型。舉個栗子，比如我知道昨天下雨，那麼今天下雨的概率就很大（狀態轉移概率）。但有時候我們並不知道當時那個地方的天氣如何，卻有可能了解到那一天雨傘的銷售量出現了提高，這樣我們就可以通過觀察雨傘的銷售量來預測天氣狀態。在這種情況中，雨傘的銷售量就是一個觀察變數，而天氣是一個隱變數，即隱藏狀態。
在股票市場中，也能遇到相似的情況。我們無法準確知道當前時刻的市場狀態，而市場狀態決定了擇時策略。但我們可以通過一系列的觀察變數來進行猜測，比如通過股票收益率、成交量、主力資金流向、融資餘額增長量等觀測數據對市場狀態（即隱藏狀態）進行猜測，並得到第二天的市場狀態預測。這樣一個基於HMM的量化模型的雛形就出來了。
三．建模步驟
對我國股市滬深300指數進行預測，假設隱藏狀態數量是6，即假設股市的狀態有6種，可以理解為「牛市上漲」、「牛市下跌」、「熊市上漲」、「熊市下跌」、「震蕩市上漲」和「震蕩市下跌」，當然實際能不能夠真正區分出這六個狀態還需要進一步驗證。
雖然我們並不知道每種狀態到底是什麼，但是通過後面的圖我們可以看出那種狀態下市場是上漲的，哪種是震蕩的，哪種是下跌的。進行預測的時候假設所有的特徵向量的狀態服從高斯分布，因此可以使用 python的hmmlearn 這個包中的 Gaussian HMM 方法進行預測了。
四．模型實現：
下面的程序均於京東量化平台實現，採用的是Python 3.5, 大家可以參考：
原文有完整代碼
閱讀原文：http://club.jr.jd.com/quant/topic/1080077
寫的再絕 人家還是覺得你用matlab的就是編程不行
google工程師的tensorflow實例代碼，theano的example…
推薦閱讀：
※作為一名數據科學家Python需要掌握到什麼程度？
※處理數據時不進行歸一化會有什麼影響？歸一化的作用是什麼？什麼時候需要歸一化？有哪些歸一化的方法？
※如何評價小米 2013 上半年業績，相對其他公司的成長有哪些可參照性？
※自己寫了一個C++矩陣計算庫，想做一下全面測試，該如何進行？