User:Trettels:Session 31 - 15APR

From Computational Statistics (CSE383M and CS395T)
Jump to navigation Jump to search

In Class

The in class activities were the same as the pre-class activity. I worked with Noah and Travis for this problem.

Main script:

mkdirdata=load('Modelselection.txt','-ascii');
t=data(:,1);
temp = data(:,2);

sm_temp=conv(temp,gausskernel(20,10),'same');
figure;
plot(sm_temp);
figure; plot(temp-sm_temp)
sd_est = temp-sm_temp;
hist(sd_est(50:end-50))
sd_est = std(sd_est(50:end-50))

plot(t,temp)
hold on; plot(sm_temp,'r'); plot(sm_temp+sd_est,'--g');plot(sm_temp - sd_est, '--g');


testfn = @(theta) (fun1(theta,temp,t,sd_est));
init = load('init.mat'); %[10 ,2*pi/450 , 10, 2*pi/500, 10, pi/(4/450),0,10,2*pi/300, 0];


options = optimset('Display','iter','MaxIter',2000, 'MaxFunEval',10000);
fn = fminsearch(testfn,init.fn,options);

[~, temp_est] = testfn(fn);
plot(t,temp_est,'k','LineWidth',1.5);

Fun1 (sum of sinusoids):

function [estlogerr, varargout] = fun1(theta,x,t,sig)

a=theta(1);
b=theta(2);
c=theta(3);
d=theta(4);
e=theta(5);
f=theta(6);
g=theta(7);
h=theta(8);
i=theta(9);
j=theta(10);

mod = (a*sin(b*t+f) + c*sin(d*t+g) + h*sin(i*t+j) + e);
estlogerr = (sum(sum(((x-mod)./sig).^2)));
if(nargout > 1)
    varargout{1} = mod;
end

Fun2 (polynomial fit):

function [estlogerr, varargout] = fun2(theta,x,t,sig)

a=theta(1);
b=theta(2);
c=theta(3);
d=theta(4);
e=theta(5);
f=theta(6);
g=theta(7);

mod = a*t.^6 + b*t.^5 + c*t.^4+ d*t.^3 + e*t.^2 + f*t.^1 + g;
estlogerr = (sum(sum(((x-mod)./sig).^2)));
if(nargout > 1)
    varargout{1} = mod;
end

The function "gausskernel" called in the main script is from the open source chronux package (www.chronux.org) and included below:

function kernel = gausskernel(R,S)
%GAUSSKERNEL       Creates a discretized N-dimensional Gaussian kernel.
%   KERNEL = GAUSSKERNEL(R,S), for scalar R and S, returns a 1-D Gaussian
%   array KERNEL with standard deviation S, discretized on [-R:R].
%   
%   If R is a D-dimensional vector and S is scalar, KERNEL will be a
%   D-dimensional isotropic Gaussian kernel with covariance matrix
%   (S^2)*eye(D), discretized on a lattice with points [-R(k):R(k)] in the
%   kth dimension.
%
%   If R and S are both D-dimensional vectors, KERNEL will be a
%   D-dimensional anisotropic Gaussian kernel on the lattice described
%   above, but with standard deviation S(k) in the kth dimension.
%
%   If R is scalar and S is a D-dimensional vector, R is treated as
%   as R*ones(D,1).
%   
%   KERNEL is always normalized to sum to 1.

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Check Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = numel(R);
D2 = numel(S);
if (((D > 1) && (~isvectord(R))) || ((D2> 1) && (~isvectord(S)))),
	error('Matrix arguments are not supported.');
end
if ((D>1)&&(D2>1)&&(D~=D2)), 
	error('R and S must have same number of elements (unless one is scalar).');
end;

if (D2>D),  D = D2;  R = R * ones(1,D); end;   % force bins/sigmas 
if (D>D2),  S = S * ones(1,D);  end;           %   to be same length,
R = R(:)';   S = S(:)';                        % and force row vectors

S(S==0) = 1e-5;  % std==0 causes problems below, 1e-5 has same effect

%%%%%%%%%%%%%%%%%%%%%%%%%%% Make the Kernel %%%%%%%%%%%%%%%%%%%%%%%%%%
RR = 2*R + 1;
for k = 1:D
	% Make the appropriate 1-D Gaussian
	grid = (-R(k):R(k))';
	gauss = exp(-grid.^2./(2*S(k).^2));  
	gauss = gauss ./ sum(gauss);

	% Then expand it against kernel-so-far ...
	if (k == 1),
		kernel = gauss;
	else	
		Dpast = ones(1,(k-1));
		expand = repmat(reshape(gauss, [Dpast RR(k)]), [RR(1:k-1) 1]);
		kernel = repmat(kernel, [Dpast RR(k)]) .* expand;
	end
end



Back to TrettelS Index