測樹學weibull擬合matlab程式，難者不會，會者不難

function [N,k,x,dbh,n_fit parmhat, parmci]=WblDbh(n,a,W)
% 根據整化的徑階，進行weibull擬合
% n : 徑階株數
% dbh : 各株的胸徑大小
% N 株數
% k 徑階數
% x 徑階胸徑
% n_fit 各徑階擬合值
% parmhat, parmci weibull擬合結果

%W = 1; %徑階距
%a = 6; %起測胸徑
dbh = [];
for i = 1:length(n)
  
    %tmp = GetDbhFrom(i,a,W)*ones(n(i),1);
    tmp = GetDbhFrom(i,0,W)*ones(n(i),1);
    dbh = [dbh;tmp];
end
N = length(dbh);
k = length(n);
%dbh = dbh - a;  % 起測胸徑平移
[parmhat parmci]= wblfit(dbh);
A=parmhat(1)
B=parmhat(2)
%figure();
f=k-3-1; %自由度
%x = W/2: W : k*W; %徑階中值
x = GetDbhFrom([1:k],0,W);
%x = GetDbhFrom([1:k],a,W);
%各徑階理論株數;
n_fit=round( N*W*B/A*power(x/A,B-1).*exp(-1*power(x/A,B)));
%n_fit = round( N*W * wblpdf(x, A,B));
figure;
plot(x +a,n,'*',x+a ,n_fit,'-o');
%chi-test
end

 

function dbh = GetDbhFrom(i,a,W)
%W = 1; %徑階距
%a = 6; %起測胸徑
%dbh 第i徑階的胸徑

dbh = W/2 + a + (i-1)*W;
end