Friday, July 1, 2011

neat matlab program I wrote

%pcabio.m

% to compute and graph PCA for stand biometrics, as well as measure distances between points!
% this file has the 10 stand metrics currently measured. The first row is
% numeric. The first column is plot names.

%------------------------------------------------------------------------
% names of files to load-- comment in what you want
%load PCAEX80.dat;
load PCABIO1984.dat
%load PCABIO2007.dat;

format short g

% creation of the file headers and column numbers. Different columns in
% early file due to no ANPP

%For 2007
% fid = fopen('PCABIO2007.dat');
% C = textscan(fid,'%s %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32');
% fclose(fid);

%for 1984
fid = fopen('PCABIO1984.dat');
C = textscan(fid,'%s %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32');
fclose(fid);

%For files that are 80
%  fid = fopen('PCAEX80.dat');
%  C = textscan(fid,'%s %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32 %f32');
%  fclose(fid);

% these are the headings that correspond to the stand metrics, in order,

%2007
%headers = {'HARDWOOD' ,'BIO' ,'PERHARD','ANPP', 'BASAL', 'NUMT', 'STEMD', 'SHANNON', 'ASSAYNO','ASSAYBIO','NUMHW','PERNUMHW','STEMDHW','PERSTEMDHW'};

%1984
headers = {'HARDWOOD' ,'BIO' ,'PERHARD','ANPP', 'BASAL', 'NUMT', 'STEMD', 'ASSAYNO','ASSAYBIO','NUMHW','PERNUMHW','STEMDHW','PERSTEMDHW','SHANNON'};

% 1980
%headers = {'HARDWOOD' ,'BIO' ,'PERHARD','BASAL', 'NUMT', 'STEMD', 'ASSAYNO','ASSAYBIO','SHANNON','NUMHW','PERHWNUM','STEMDHW','PERSTEMDHW'};

% create a matrix that does not include the first row of PCAEX07 because it
% is text for the plot numbers (see above string in %s)

%X = PCAEX80(:,2:14);
X = PCABIO1984(:,2:15);
%X = PCABIO2007(:,2:15);

%------------------------------------------CONDUCTING PCA!!!------------

% if there are no NaN's!

[PC, scores, latent] = princomp(zscore(X));

% if NaN's are included, you must centralize and standardize each row
% individually, so comment in this code!

% % subtract out the mean
% for i = 1:length(X)
% thisRow = X(i,:);
% thisRow(isnan(thisRow)) = [];
% stdX(i) = std(thisRow);
% meanX(i) = mean(thisRow);
% centerX = X - repmat(meanX,10,1);
% centerStandardX = centerX./repmat(stdX,10,1);
% end


% % conduct principle components analysis
% [PC, scores, latent] = princomp(centerStandardX);

%----------------DISPLAYING THE CUMULATIVE OUTPUT OF THE VARIANCE!---

account = cumsum(latent)./sum(latent);
disp(account);

% name eigenvectors to the vectors that you are about.
vect1 = scores(:,1);
vect2 = scores(:,2);

% plot the principle components (FACTOR ANALYSIS PLOT) with the plot
% numbers-- must make sure these are STRINGS and not FLOATS/INTS
figure(1)
plot(vect1,vect2,'b.')
xlabel('1st Principal Component')
ylabel('2nd Principal Component')
title('Principal Components of Stand Biometerics')
text(vect1+0.1, vect2, C{1}, 'FontSize',8);

% plot a scree plot of cumulative variance explained
figure(2)
pareto(latent)
xlabel('Principal Component')
ylabel('Variance Explained (%)')
title('Scree plot of variance explained and CDF')

% plot a biplot for the first two axes... can be done for 3 as well, if
% needed, but is hard to interpret in my opinion. if disp(account)'s second
% entry is bigger than about 0.55, I don't do it.

figure(3)
biplot(PC(:,1:2), 'scores',scores(:,1:2),'varlabels',headers')


% for finding the distance between the sets of points and printing them out
% by labels! Will be used later in geostats

b = zeros(length(vect1),2);
for i = 1: length(vect1)
    b(i,1) = vect1(i);
    b(i,2) = vect2(i);
end

k = 0;
for i = 1:length(b)
    for j = 1:length(b)
        if j~=i
            k=k+1;
            index(k,:) = [i j];
        end;
    end;
end;

dist = ((b(index(:,2),1)-b(index(:,1),1)).^2 + (b(index(:,2),2)-b(index(:,2),2)).^2).^0.5;

% if you want to print it to the scree-- WARNING IT WILL OVERLOAD YOUR
% MEMORY IF YOU ARE RUNNING MUCH ELSE!

% fprintf(' #1 #2 distance\n');
% fprintf('%3d %3d %10.6f\n', ...
%    [index(:,1),index(:,2),dist].');
% fprintf('Total %11.6f\n',sum(dist));

m = [index(:,1),index(:,2),dist];

% output the distance matrix information to a file that saves in the Matlab
% directory. I think this format will be adaptable well to SPSS, but it
% would be better for Matlab to have in matrix form!

%dlmwrite('biometric_distances_80.csv',m,',');
dlmwrite('biometric_distances_804.csv',m,',');
%dlmwrite('biometric_distances_07.csv',m,',');

Sunday, June 26, 2011

This is kind of an elegant little solution

This code can be used to find the distance between x locations (stored in vect1) and y locations (stored in vect2). I modified this from something I found online. Works pretty good and I like the out put, also optimal in terms of speed.
Using this to construct similarity ranks! win!

% for finding the distance between the sets of points and printing them out
% by labels!

b = zeros(length(vect1),2);
for i = 1: length(vect1)
    b(i,1) = vect1(i);
    b(i,2) = vect2(i);
end

k = 0;
for i = 1:length(b)
    for j = 1:length(b)
        if j~=i
            k=k+1;
            index(k,:) = [i j];
        end;
    end;
end;

dist = ((b(index(:,2),1)-b(index(:,1),1)).^2 + (b(index(:,2),2)-b(index(:,2),2)).^2).^0.5;
m = [index(:,1),index(:,2),dist];

Saturday, June 25, 2011

Soil depth

I expected this but it still seems funny that soils are deepest at low... and high... elevations. 

It's because on this landscape the ridges are relatively flat, and are located between precipices out of the extent of the study area that are eroding soils onto them. Still. Funny. 

Friday, June 24, 2011

best help manual I've ever seen

http://www.box.net/shared/mugpnnjy21

Ikea couldn't have done it better.

How to read text for matlab

fid = fopen('biomatlabels.dat');
C = textscan(fid,'%s %f32 %f32 %f32 %f32 %f32 %f32 %f32');
fclose(fid);

The first line opens the file. You must save it as a .dat file!

The second line scans the file. it is also important to name it fid.
%s means "string" (like a word)
%f32 means "number"(in the broadest sense-- there are options for decimals, too)

The third line closes the file. When you want to call up one of your labels, call it as
C{1} or C{2}-- it saves as a cell array so brackets are key!

Wednesday, June 8, 2011

a code for stepwise regression with twelve possible fitting parameters

% litter.m


file = importdata('biomass_workbook2.csv',',');


% names of information | bio in 1980 | percenthardwood in 1980 |
% hardwood bio 1980|...etc.
bio80 = file(:,1);
perhard80 = file(:,2);
hardbio07 = file(:,3);
perhard07 = file (:,4);
litterdat = file(:,5);
anpp1 = file(:,7);
basal07 = file(:,8);
anpp6 = file(:,9);
bio07 = file(:,11);
num_tree = file(:, 12);
stemden = file(:,13);
logherb = file(:,14);
hwbio80 = file (:,15);


% poly fit finds the parmaters (pi) and the norms of the residuals (si) for
% linear and poly combinations of variables


[p1,s1] = polyfit(bio80, litterdat, 1);
[p2,s2] = polyfit(perhard80, litterdat, 1);
[p3,s3] = polyfit(hardbio07, litterdat, 1);
[p4,s4] = polyfit(perhard07, litterdat, 1);
[p5,s5] = polyfit(anpp1, litterdat, 1);
[p6,s6] = polyfit(basal07, litterdat, 1);
[p7,s7] = polyfit(anpp6, litterdat, 1);
[p8,s8] = polyfit(bio07, litterdat, 1);
[p9,s9] = polyfit(num_tree, litterdat, 1);
[p10,s10] = polyfit(stemden, litterdat, 1);
[p11,s11] = polyfit(logherb, litterdat, 1);
[p12,s12] = polyfit(hwbio80, litterdat, 1);


% polyval finds the correlation matrix
output1 = polyval(p1, bio80);
corr1 = corrcoef(litterdat,output1);


output2 = polyval(p2, perhard80);
corr2 = corrcoef(litterdat,output2);


output3 = polyval(p3, hardbio07);
corr3 = corrcoef(litterdat,output3);


output4 = polyval(p4, perhard07);
corr4 = corrcoef(litterdat,output4);


output5 = polyval(p5, anpp1);
corr5 = corrcoef(litterdat,output5);


output6 = polyval(p6, basal07);
corr6 = corrcoef(litterdat,output6);


output7 = polyval(p7, anpp6);
corr7 = corrcoef(litterdat,output7);


output8 = polyval(p8, bio07);
corr8 = corrcoef(litterdat,output8);


output9 = polyval(p9, num_tree);
corr9 = corrcoef(litterdat,output9);


output10 = polyval(p10, stemden);
corr10 = corrcoef(litterdat,output10);


output11 = polyval(p11, logherb);
corr11 = corrcoef(litterdat,output11);


output12 = polyval(p2, hwbio80);
corr12 = corrcoef(litterdat,output12);


% M is a matrix to hold thePearson's correlation coefficient
m = zeros(12,1);
m(1) = corr1(1,2);
m(2) = corr2 (1,2);
m(3) = corr3 (1,2);
m(4) = corr4 (1,2);
m(5) = corr5 (1,2);
m(6) = corr6 (1,2);
m(7) = corr7 (1,2);
m(8) = corr8 (1,2);
m(9) = corr9 (1,2);
m(10) = corr10 (1,2);
m(11) = corr11 (1,2);
m(12) = corr12 (1,2);


% yh(1) holds the fit of the predicted 


yh1 = polyval(p1,bio80);
yh2 = polyval(p2, perhard80);
yh3 = polyval(p3, hardbio07);
yh4 = polyval(p4, perhard07);
yh5 = polyval(p5, anpp1);
yh6 = polyval(p6, basal07);
yh7 = polyval(p7, anpp6);
yh8 = polyval(p8, bio07);
yh9 = polyval(p9, num_tree);
yh10 = polyval(p10, stemden);
yh11 = polyval(p11, logherb);
yh12 = polyval(p12, hwbio80);


% resid is what is not predicted


resid1 = litterdat - yh1;
resid2 = litterdat - yh2;
resid3 = litterdat - yh3;
resid4 = litterdat - yh4;
resid5 = litterdat - yh5;
resid6 = litterdat - yh6;
resid7 = litterdat -yh7;
resid8 = litterdat - yh8;
resid9 = litterdat - yh9;
resid10 = litterdat - yh10;
resid11 = litterdat - yh11;
resid12 = litterdat - yh12;


rss1 = sum((yh1- mean(litterdat)).^2);
tss1 = sum((litterdat - mean(litterdat)).^2);


% Rsq asks how much better predictedis thanmean
Rsq1 = 1-sum(resid1.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq2 = 1-sum(resid2.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq3 = 1-sum(resid3.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq4 = 1-sum(resid4.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq5 = 1-sum(resid5.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq6 = 1-sum(resid6.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq7 = 1-sum(resid7.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq8 = 1-sum(resid8.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq9 = 1-sum(resid9.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq10 = 1-sum(resid10.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq11 = 1-sum(resid11.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq12 = 1-sum(resid12.^2)/sum((litterdat-mean(litterdat)).^2);


% n is a matrix of R sq
n= zeros(12,1);
n(1) = Rsq1;
n(2) = Rsq2;
n(3) = Rsq3;
n(4) = Rsq4;
n(5) = Rsq5;
n(6) = Rsq6;
n(7) = Rsq7;
n(8) = Rsq8;
n(9) = Rsq9;
n(10) = Rsq10;
n(11) = Rsq11;
n(12) = Rsq12;


% mat is a matrix of predictors
mat = zeros(length(litterdat),12);
mat(:,1:4) = file(:,1:4);
mat(:,5:7) = file(:,7:9);
mat(:,8:12) = file(:,11:15);


% this will do a stepwise fit. b isparameters, se is residuals, pval is
% values.. in model shows the process of gettingthere, statsis stats,next
% step and history are more dtailed
[b, se, pval, inmodel, stats, nextstep, history] = stepwisefit(mat, litterdat);


% model_fit is parameters
model_fit = inmodel.*b';
model_fit = model_fit';


% model 2is the predicted
model2 = zeros(length(mat),1);
model2= 4.4443 + model_fit(4).*mat(:,4) + model_fit(7)*mat(:,7) + model_fit(8).*mat(:,8);


% calculateresiduals
resid_step = litterdat-model2;


% calculate r2
Rsq_step = 1-sum(resid_step.^2)/sum((litterdat-mean(litterdat)).^2);


% calculate a fitting "line" to plot
[p_step, s_step] = polyfit(model2, litterdat,1);
a = p_step(1);
b = p_step(2);
Vcalc = a.*model2 + b;


% plot
figure(1)
plot(model2,litterdat,'k.');
hold on
plot(model2, Vcalc,'r-');
xlabel('Stepwise litterfall model')
ylabel('Litterfall data')
title ('Fitness testing for litterfall model');

Thursday, June 2, 2011

Processing Tutorials

 I'm working on my data visualization skills here a bit, and I thought I'd start trying to figure out how to work with Processing. Although ultimately I think VELMA will switch languages, Processing does have the advantage of NOT being in C and looking nice.
So to the right is the start of an example from
Data Visualization. Since I'm a total n00b I thought, I wonder what happens when we change the size of the image.

As you can see here, changing the first number cuts it off from the right to left, and the second number from the bottom to the top.

If we change the image numbers, we see that going from 0 to 100 shifts the image within the background.
This, I suppose, is like being in "data mode" in Arc.