Wed 19 Nov 2014 08:42:33 AM UTC, comment #2:
I have no idea why the rest of my post was truncated. Here it is again:
"##S mapping
pkg load statistics
cd C:\Octave\Octave-3.8.2
load 'logistic.dat' #load up and pre-process a timeseries of data from the
logisic map to use as an example
#as it is chaotic.
data=flipud(logistic); #flip so that most recent at top, so not going
backwards.
#load 'tinkerbell.dat';
#data=flipud(tinkerbell);
#may need to then one minus the data if doing
#one step projection, or Tp minus if doing Tp steps.
#define m and d, theta and Tp (last 3 will be in function when done)
#d=2;
#Tp=1
#theta=2
#latest = 90
#earliest = 150
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;
Yhat = (A\B).*Xt;
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
"
I then pasted in the following code, which does have an error in it (it
doesn't make sense for i to start at 0), and this led to the whole thing
crashing. (It doesn't crash if changed to be i=1:10).
"
clvector=1
for i=0:10
[yt,cl] = smap(data,150,90,i,1,4);
clvector = [clvector;cl];
endfor
clvector = clvector(2:(length(clvector)))
axis = 0:(length(clvector)-1)
plot (axis, clvector, 'x')
"
There you go, that's it. Logistic.dat is the data file used in the code, it's just a time series.
|