

datadir='C:\Users\klidke\Documents\Teaching\480BiophysicsFall2010\MATLAB\';
data=load([datadir 'SimDiff.mat']);

x=data.x;
y=data.y;

T=length(x);
pixelsize=.1;% micron
timestep=.01;% second

%I'm going to do everything in pixels and frames, then convert at the end
%to physical units.

%Make MSD
%I'm not re-using points for a specific deltat
Tmax=round(T*.25) %make MSD with deltat only up to 25% of total time
MSD=zeros(Tmax,1);
N=zeros(Tmax,1);
deltat=1;
while deltat<=Tmax %this loops over deltat
    tmp=0;
    ii=1;
    while (ii+deltat)<=T % This loops over the start point and stops when past the end
        N(deltat)=N(deltat)+1; %counts how many jumps are used
        r2=(x(ii+deltat)-x(ii))^2+(y(ii+deltat)-y(ii))^2;
        tmp=tmp+r2;
        ii=ii+deltat;
    end
    MSD(deltat)=tmp/N(deltat); 
    deltat=deltat+1;
end

%adjust for physical units
MSD_p=MSD*pixelsize^2;
deltat_p=(1:Tmax)'*timestep;

figure
plot(deltat_p,MSD_p)
xlabel(' \Delta t (s)')
ylabel('MSD(\Delta t)')

%% Fitting the MSD

% There are many ways to perform the fit.  I'm using an approach that does
% not require any toolboxes. 

%weighting:
%In a chi-squared minimization, each point is weighted with 1/variance. In our case, the
%variance is proportional to number of points, N,  used for each MSD bin.
%The error is proportional to D*deltat (MSD=4Dt), so variance is proportional to 
% (D*deltat)^2. D is a fit parameter, and is the same at each MSD point,
% therefore, we just need the deltat^2.  
%ref: Qian, Biophyiscal Journal, 1991.   

V=N.*deltat_p.^2;

%I'm using fminsearch with an inline function
%this is sum( (data-model)^2/Variance )
%model is a+b*deltat_p+c*deltat_p^2
%since fminsearch takes a parameter vector, X=[a b c];

LSQMSDfit=inline('sum( (MSD-( X(1)+X(2)*deltat_p+X(3)*deltat_p.^2)).^2./V)','X','MSD','deltat_p','V');
options=optimset;
X0=[.0001 .01 .005]; %these are the starting guesses
% X0 and the inputs after 'options' get passed to the function 'func'
[out]=fminsearch(LSQMSDfit,X0,options,MSD_p,deltat_p,V);

%check fit
%out=X0 %for checking starting guess
model=out(1)+out(2)*deltat_p+out(3)*deltat_p.^2;
figure
plot(deltat_p,MSD_p);
hold on
plot(deltat_p,model,'r');
xlabel(' \Delta t (s)')
ylabel('MSD(\Delta t)')

Pos_err=sqrt(out(1)/2)  %microns
D_found=out(2)/4        %micron^2/s
V_found=sqrt(out(3))    %microns/s

% The values that generated the simulated data were:
%Perr   =    0.003000000000000 microns
%v      =    0.086023252670426 microns/s
%D      =    0.001500000000000 microns^2/s












