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!

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');

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.

so much more to know

to be a great modeler, I must know more math.
to be a decent ecologist, I must know more ecology.
to be a great spatial analyst, I must become better at programming.
to understand forests better, I must work outside.
to be smart like everyone else, I must know hydrology, fire ecology, landscape ecology, statistics, bioenergetics, physics, chemistry, ecophysiology...

in short, I am behind. it's daunting. I know I need a year to just nose to the grind and learn. I don't know if I can afford it. But I think I need it. I think I could come out with some mediocre stuff in the meanwhile to keep me afloat-- you know, just helping out with spatial stuff here and there, valuation, etc.


One of many things I need to learn is Bayes theorem. I don't get it at all. I've never learned it in a class or talked to anyone about it. I've tried to read about it but it's light years ahead of me. I have it memorized, but it doesn't "click" without practice. This one passage about it, though, I found pretty helpful, so I thought I'd share.

From an article in Ecology by Subhash R. Lele (2010, vol. 91 (12)) - Big Fancy Models:


On the other hand, if the posterior distribution converges to a nondegenerate distribution, it implies non-estimability of the parameters. This nondegenerate distribution can be, and usually is, different than the prior distribution; there can be ‘‘Bayesian learning’’ without identifiability.
Consider a simple example. Let Yi conditional Mu ~ N (Mu, sigma^2)
and let Mu i ~ N (l, tau^2)
Then it is obvious that Yi ~N (l, sigma^2 + tau^2). The parameters sigma^2 and tau^2 are individually nonidentifiable. Suppose we put priors sigma^2~Unif(0, 100) and tau^2~ Unif(0, 100). Suppose the truth is such that sigma^2 + tau^2=  10. Then the marginal posterior distributions for sigma^2 and tau^2 necessarily get concentrated on the interval (0, 10) as the sample size increases. Their joint distribution will be concentrated along the diagonal of the square defined by the coordinates (0, 0), (0, 10), (10, 10), and (10,0). This distribution is different than the prior distribution. Thus, there is ‘‘Bayesian learning’’ but clearly existence of Bayesian learning does not imply that the parameters are identifiable or even that legitimate inferences can be drawn about the parameters for which Bayesian learning happens. If a part of the model is non-identifiable, it can make estimators of other parameters inconsistent. They converge to a single, but wrong point...



Ecologists know a great deal about the processes. While constructing mathematical models, they have a strong and admirable desire to include all the nuances. Unfortunately the data are not always informative enough to conduct inferences on all the complexities of the model. As a consequence, either the model parameters become non-identifiable or non-estimable. If estimation is possible, estimates tend to be extremely uncertain with large standard errors, thus precluding their use in effective decision making. I would urge ecologists to establish identifiability of the parameters in their models before conducting any scientific inferences...

A program for calculating BAHA, Number of Trees, Biomass

This program can be used to calculate basal area per hectare, number of trees, and biomass/ha from an input file. It is modified from a program I wrote (With help!) this summer

% live_doug.csv


%---THE REAL PROGRAM BEGINS HERE!!!----------------------------------------%


% Import the data from donfile.csv into the variable "bio." This is the biomass
% data calculated directly with the equations from the forest inventories
% in the field. 
% | YEAR | TREE NO | DBH | BIO | TRANS | PLOT | PLOTID
% the domains for these inputs are years, index, index, index,
% concat(indexindex), concat (yearsindex),Mg, Mg, Mg, Mg, and Mg,
% respectively. Each row represents and individual tree.


FILE = 'live_doug_2007.csv';


% Import the LIDAR data which is courtesy of Keith. This data is in the
% following format:
%   | TRANSECT  | PLOT | PLOTID |AREA OF ELLIPSE 
% the domains for these inputs are index, index, concat(indexindex), and Ha, 
% respectively. The slope, aspect, and Ha of each ellipse was calculated
% using ArcGIS-- it is not part of this program. 


FILE_2 = 'ws01_attrib_lidar.csv';




%giving the files some easier names to use in the program
bio = importdata(FILE,',');
hectares = importdata(FILE_2,',');


% Find the unique values in the concatenated years & plot number
% column, and store their positions in upos
[uval,upos] = unique(bio(:,7));


% Define an Nx8 matrix to hold the sums down columns
% for every unique position
sum_bio = zeros(length(upos),6);


% Copy the years as column 1 of the sum matrix, the plot # as col 2,
% and the unique ID as column 3 of the sum matrix
sum_bio(:,1) = bio(upos,1);
sum_bio(:,2) = bio(upos,7);


% Sum down the columns of bio according to the elements of upos
% then divide by the appropriate element of the hectares matrix,
% and store the results in the sum_bio matrix


% sum_bio looks like (all sums divided by hectares):
% year-1966 | plot | ID | sum 1 | sum 2 | sum 3 | sum 4 | sum total
sum_bio(1,3:4) = sum(bio(1:upos(1),3:4))./hectares(find(hectares(:,3) == sum_bio(1,2)),4);
for j = 1:(length(upos)-1);
    sum_bio(j+1,3:4) = sum(bio((upos(j)+1):upos(j+1),3:4))./hectares(find(hectares(:,3) == sum_bio(j+1,2)),4);
end


% number of trees
sum_bio(1,5) = upos(1);


for i = 2:length(upos-1)
sum_bio(i,5) = upos(i)-upos(i-1);
end


%BAHA


sum_bio(1,6) = sum_bio(1,3)*0.00007854;


for i = 2:length(upos-1)
sum_bio(i,6) = sum_bio(i,3).*0.00007854;
end


% export all information to a separate CSV file in MatLab directory
%dlmwrite('bioassay2007.csv', sum_bio, ',');

when math doesn't speak english, part 2

I found this little math gem in a paper I"m reading about Canonical Correspondence Analysis, a technique which has the potential to be a nice "balance" between going completely geostatistical and sticking to the ecologist-friendly realm of P-Value-Ville.

So this little joy here simply means "Mean"-- it's the definition of the mean. Does that look like a "mu" to you? I have to pretend sometimes. 

The paper this is form is 

Canonical Correlation Analysis: An overview with application to learning methods. David R Hardoon, S. Szedmak, and John Shawe-Taylor. Technical Report: CSD-TR-03-02. 2003.