function [Yhat, closeness, weights, FAR] = smap(data,earliest,latest,d,Tp,theta) x=data(latest:earliest,:); m=length(x); #always rescale so Xt(0) = 1 //try not doing this #x = x.+(1-(x(m,:))); #embed the time series: m = length(x); embedded = ones((m-d),1); #presize a matrix for i = 1:d col = x([i:(m-(d-(i-1)))],1); embedded = [embedded, col]; endfor # this embedded has the original column of 1's, so chop them off: # Actually, don't! embedded= embedded(:,2:(d+1)); # embedded is now a matrix of d+1 columns and n-d-1 rows. [n,dim] = size(embedded); # get size of the embedded time series matrix #calculate distances between all points to get distbar #can do this using pdist (statistics package) so need to load statistics if not already done distbar = mean(pdist(embedded)) #now need the distance between latest simplex and others: #Assumes most recent at the top Xt = embedded(1,:); #Need Yi: Yi = embedded(2:Tp:n-Tp,:); #And Xi: Xi = embedded(Tp+2:Tp:n,:); #FAR is vector of distances between latest known point and every other point in library #i.e Work out (||Xi-Xt||) far = Xi.-Xt; FAR = sqrt(sum(far.^2,2)); #calculate weights These are not normalized, so divide by sum weights weights = exp(-theta.*FAR./distbar); weights = weights./(sum(weights)); #calculate Yhat: yhat = A inverse B Xt #where A = weight(||xi-Xt||)Xi and B = weight(||Xi-Xt||)Yi #Xi is each point in the library, Xt is the predictee (point) and Yi is where #Xi ended up at. #Yhat = mean(weights.*Yi); #Yhat = Yhat(1,1) B = weights.*Yi; A = weights.*Xi; if (~any(isnan(A)) && ~any(isnan(B))) Yhat = (A\B).*Xt; else warning("There are unexptected values in A or B") Yhat = Xt; endif Yhat = sum(Yhat); Yhat = Yhat(1,2) #now get the actual true value for the point from the data: tru = data(1:(latest-1),:); #always rescale so Xt(0) = 1 //try not doing this # tru = tru.+(1-(x(m,:))); #Difference from actual point Yt: Yt = tru(length(tru)); closeness = sqrt((Yhat.-Yt)^2); endfunction