In [186]:
time = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
pos = [1, 1.1, 0.95, 0.9, 0.2, 0.1, 0, 0.1, 0.01, 0.2,0.1];
figsize(6,4)
plot(time,pos); title('Dummy data');
plt.xlabel('time'); plt.ylabel('signal');
Let's define the biasing parameters
In [187]:
So = 0; nu = 0
In [188]:
firstDwellMean = mean(pos);
plot(time,pos);
plot([time[0],time[10]],[firstDwellMean,firstDwellMean]);
title('Single dwell fit'); xlabel('time'); ylabel('signal');
In [189]:
print 'Sigma squared is: ' + str(var(pos))
print 'Standard deviation is: ' + str(sqrt(var(pos)))
firstDwellSigmaSquared = var(pos);
In the figure below, the blue line shows the raw data, the green line the mean of the data, and the red lines the standard deviation range.
In [190]:
plot(time,pos);
plot([time[0],time[10]],[firstDwellMean,firstDwellMean]);
plot([time[0],time[10]],[[firstDwellMean+sqrt(var(pos))],[firstDwellMean+sqrt(var(pos))]], 'r');
plot([time[0],time[10]],[[firstDwellMean-sqrt(var(pos))],[firstDwellMean-sqrt(var(pos))]], 'r');
title('Single dwell fit'); xlabel('time'); ylabel('signal');
Let's calculate the SICP for this fit.
In [191]:
firstDwellSICP = -1*(log(2*pi)+1)+log(11)+(11+nu-1-1)*log(11*firstDwellSigmaSquared+So)-(11+nu-1-3)*log(11+nu-1-3)-log(11+nu-1-3)
print "The first dwell's SICP is: " + str(firstDwellSICP)
In [279]:
print "The SICP for each possible step location is: "
singleStepSICP=[] # this is a list of every SICP for each single step possibility
singleStepMean=[] # this is a list of every set of dwell means for each single step possibility, has form [[],[],...]
singleStepSigmaSquared=[] # this is a list of every set of dwell sigma squared for each single step possibilty, has form [[],[],...]
for stepLocation in range(1, len(pos)): # iterate through each possible single step location
firstDwellMean = mean(pos[0:stepLocation])
secondDwellMean = mean(pos[stepLocation:len(pos)])
singleStepMean.append([firstDwellMean,secondDwellMean])
firstDwellSigmaSquared = var(pos[0:stepLocation])
secondDwellSigmaSquared = var(pos[stepLocation:len(pos)])
singleStepSigmaSquared.append([firstDwellSigmaSquared,secondDwellSigmaSquared])
singleStepOverallSigmaSquared = (len(pos[0:stepLocation])*firstDwellSigmaSquared+len(pos[stepLocation:len(pos)])*secondDwellSigmaSquared)/len(pos)
currentSICP = -2*(log(2*pi)+1)+log(firstDwellMean)+log(secondDwellMean)+(len(pos)+nu-2-1)*log(len(pos)*singleStepOverallSigmaSquared+So)-(len(pos)+nu-2-3)*log(len(pos)+nu-2-3)-log(len(pos)+nu-2-3)
singleStepSICP.append(currentSICP)
print str(singleStepSICP)
print "So the optimal step location for the first step is at index " + str(singleStepSICP.index(min(singleStepSICP))) + " corresponding to SICP value " + str(min(singleStepSICP))
We can extract the optimal dwell means for the single step calculation (located at index 3) very easily
In [280]:
print "The respective dwell means are: " + str(singleStepMean[3])
singleStepFit = []
for dataPoint in range(0,singleStepSICP.index(min(singleStepSICP))+1):
singleStepFit.append([time[dataPoint],singleStepMean[3][0]])
for dataPoint in range(singleStepSICP.index(min(singleStepSICP)),len(pos)):
singleStepFit.append([time[dataPoint],singleStepMean[3][1]])
singleStepFit
plot(time,pos)
plot([x[0] for x in singleStepFit],[x[1] for x in singleStepFit]);
The plot above shows the optimal single step fit as determined by the SICP. We need to compare this to the optimal no step fit to make sure that our added step fits better than the previous fit
In [294]:
if min(singleStepSICP) < firstDwellSICP:
print "The single step fits better than having no step " + "because it has a lower SICP"
else:
print "You should never see this"
In [ ]:
print "The SICP for each possible additional step location is: "
doubleStepSICP=[] # this is a list of every SICP for each two step possibility
doubleStepMean=[] # this is a list of every set of dwell means for each two step possibility, has form [[],[],...]
doubleStepSigmaSquared=[] # this is a list of every set of dwell sigma squared for each single step possibilty, has form [[],[],...]
for stepLocation in range(1, len(pos)): # iterate through each possible single step location
if stepLocation not in singleStepSICP.index(min(singleStepSICP)):
firstDwellMean = mean(pos[0:stepLocation])
secondDwellMean = mean(pos[stepLocation:len(pos)])
singleStepMean.append([firstDwellMean,secondDwellMean])
firstDwellSigmaSquared = var(pos[0:stepLocation])
secondDwellSigmaSquared = var(pos[stepLocation:len(pos)])
singleStepSigmaSquared.append([firstDwellSigmaSquared,secondDwellSigmaSquared])
singleStepOverallSigmaSquared = (len(pos[0:stepLocation])*firstDwellSigmaSquared+len(pos[stepLocation:len(pos)])*secondDwellSigmaSquared)/len(pos)
currentSICP = -2*(log(2*pi)+1)+log(firstDwellMean)+log(secondDwellMean)+(len(pos)+nu-2-1)*log(len(pos)*singleStepOverallSigmaSquared+So)-(len(pos)+nu-2-3)*log(len(pos)+nu-2-3)-log(len(pos)+nu-2-3)
singleStepSICP.append(currentSICP)
print str(singleStepSICP)
print "So the optimal step location for the first step is at index " + str(singleStepSICP.index(min(singleStepSICP))) + " corresponding to SICP value " + str(min(singleStepSICP))
In [290]:
range(0,len(pos))
Out[290]:
A fit is fully specified by:
In [289]:
"fit_"+str(3)+"steps" = 3
With non-stationary noise we need to slice the dataset into small enough chunks such that the data w/in each chunk can be assumed stationary. A typical dataset will have time, force, and position. The noise will change as a function of the force. Slice the dataset by force.
In [ ]: