maxminest.m

低矮房屋风压系数、风荷载计算分析matlab程序

function [max_est min_est] = maxminest( record, duration_ratio )

% The function, written by Joseph A. Main is a modified version of a 
% function previously written by Fahim Sadek, based on the procedure 
% described in the following paper:
% 
% Sadek, F. and Simiu, E. (2002). "Peak non-gaussian wind effects for
% database-assisted low-rise building design." Journal of Engineering
% Mechanics, 128(5), 530-539.
%
% The function computes expected maximum and minimum values of time series
% by estimating probability distributions for the peaks.

% Each row of 'record' is a time series.
% 'duration_ratio' = (duration for expected peaks)/(duration of record)

if nargin==1 || isempty( duration_ratio )
    duration_ratio = 1;
end

plot_on = 0; % turns plotting on (1) [for diagnostics] or off (0)


n_cdf_pk = 1000;     % number of points in CDF for peaks
cdf_pk_min = 0.0005; % minimum value of CDF for peaks
cdf_pk_max = 0.9995; % maximum value of CDF for peaks
cdf_pk = linspace(cdf_pk_min,cdf_pk_max, n_cdf_pk); % linearly spaced CDF values for peaks

rsize = size(record);

max_est = zeros(rsize(1),1);
min_est = zeros(rsize(1),1);

for i=1:rsize(1)
    x = record(i,:);
    n = length(x);

    mean_x = mean( x );
    std_x = std(x, 1); % biased estimate
    skew_x = sum((x - mean_x).^3)/ (n*std_x^3); % biased estimate (Matlab default for 'skewness')
    
    x = x*sign(skew_x); % change sign (if necessary) to have positive skewness
    
    sort_x = sort( x );
    CDF_x = (1:n)/(n+1);
    
    % Decide where to split the CDF (the lower portion will be fit with a
    % normal distribution and the upper part with a gamma distribution):
    CDF_split_list = linspace(0.2,0.5,10); % search from the 20th percentile to the 50th percentile of the values
    x_split_list = interp1(CDF_x, sort_x, CDF_split_list);
    norm_PPCC_list = zeros(size(x_split_list));
    
    s_norm = stdnorminv(CDF_x); % standard normal variate
    
    for j = length(x_split_list):-1:1
        
        % Obtain the Normal Distribution Parameters for lower portion of CDF
        ind_lo_j = find(sort_x<x_split_list(j));
        x_lo_j = sort_x(ind_lo_j);
        n_lo_j = length(x_lo_j);
        CDF_lo_j = CDF_x(ind_lo_j);
        s_norm_lo_j = s_norm(ind_lo_j); % standard normal variate
        mean_s_norm_lo_j = mean(s_norm_lo_j);
        mean_x_lo_j = mean(x_lo_j);

        sigma_lo_list(j) = (sum(s_norm_lo_j(:).*x_lo_j(:))-n_lo_j*mean_s_norm_lo_j*mean_x_lo_j)/(sum(s_norm_lo_j.^2)-n_lo_j*mean_s_norm_lo_j^2);
        mu_lo_list(j) = mean_x_lo_j - sigma_lo_list(j)*mean_s_norm_lo_j;
        x_lo_fit_j = mu_lo_list(j) + sigma_lo_list(j)*s_norm_lo_j;
        
        norm_PPCC_list(j) = sigma_lo_list(j)*std(s_norm_lo_j)/std(x_lo_j);

        if plot_on
            figure(1)
            plot(s_norm_lo_j,x_lo_j,'.',s_norm_lo_j,x_lo_fit_j,'-');
            title(['x_split: ' num2str(x_split_list(j)) '    PPCC: ' num2str(norm_PPCC_list(j))]);
            pause;
        end

        % Split the CDF at the value that maximizes the Probability Plot
        % Correlation Coefficient:
        if norm_PPCC_list(j)==max(norm_PPCC_list)
            norm_PPCC = norm_PPCC_list(j);
            x_split = x_split_list(j); 
            x_lo = x_lo_j;
            CDF_lo = CDF_lo_j;
            sigma_lo = sigma_lo_list(j);
            mu_lo = mu_lo_list(j);
            s_norm_lo = s_norm_lo_j;
            x_lo_fit = x_lo_fit_j;
        end
    end
    if plot_on
        figure(1)
        plot(x_split_list,norm_PPCC_list,'.',x_split,norm_PPCC,'o')
        pause;
    end
    
    % Estimate gamma distribution parameters for upper portion of CDF:
    ind_hi = find(sort_x>=x_split);
    x_hi = sort_x(ind_hi);
    n_hi = length(x_hi);
    CDF_hi = CDF_x(ind_hi);
    
    mean_x_hi = mean(x_hi);
    std_x_hi = std(x_hi);

%    gamma_MOM = (2.0/skew_x)^2; % estimate of gamma from Method of Moments (for entire time series, not just upper portion)
    gamma_min = 1;
    gamma_max = 125;
    n_gamma = 15;
    gamma_list = logspace(log10(gamma_min),log10(gamma_max),n_gamma);
    gam_PPCC_list = zeros(size(gamma_list));
%     if gamma_MOM<gamma_max
%         gamma_list = union(gamma_list,gamma_MOM);
%     end
    for j = length(gamma_list):-1:1
        % Obtain the Normal Distribution Parameters for lower portion of CDF
        s_gam_j = stdgaminv(CDF_hi(:), gamma_list(j)); % standard variate
        mean_s_gam_j = mean(s_gam_j);

        beta_hi_list(j) = (sum(s_gam_j(:).*x_hi(:))-n_hi*mean_s_gam_j*mean_x_hi)/(sum(s_gam_j.^2)-n_hi*mean_s_gam_j^2);
        mu_hi_list(j) = mean_x_hi - beta_hi_list(j)*mean_s_gam_j;

        gam_PPCC_list(j) = beta_hi_list(j)*std(s_gam_j)/std(x_hi);
        x_hi_fit_j = mu_hi_list(j) + beta_hi_list(j)*s_gam_j;
        if gam_PPCC_list(j)==max(gam_PPCC_list)
            gam_PPCC = gam_PPCC_list(j);
            gamma_hi = gamma_list(j);
            beta_hi = beta_hi_list(j);
            mu_hi = mu_hi_list(j);
            s_gam = s_gam_j;
            x_hi_fit = x_hi_fit_j;
        else
            break; % stop searching once the PPCC starts to decrease
        end
        if plot_on
            figure(2)
            plot(s_gam_j,x_hi,'.',s_gam_j,x_hi_fit_j,'-');
            title(['gamma: ' num2str(gamma_list(j)) '    PPCC: ' num2str(gam_PPCC_list(j))]);
            pause;
        end
    end
    if plot_on
        figure(2)
        plot(gamma_list(find(gam_PPCC_list)),gam_PPCC_list(find(gam_PPCC_list)),'.',gamma_hi,gam_PPCC,'o')
        pause;
        figure(3);
        subplot(2,1,1)
        plot(sort_x,CDF_x,'.',x_hi_fit,CDF_hi,'-',x_lo_fit,CDF_lo,'-');
        subplot(2,2,3)
        plot(s_norm_lo,x_lo,'.',s_norm_lo,x_lo_fit,'-');
        subplot(2,2,4)
        plot(s_gam,x_hi,'.',s_gam,x_hi_fit,'-');
    end
   
    % Compute the mean zero upcrossing rate of process y(t) with standardized
    % normal probability distribution.
    x_u = sort_x( find(abs(CDF_x - 0.5)==min(abs(CDF_x - 0.5))) ); % upcrossing level
    Nupcross = length(find( x(2:end)>=x_u & x(1:end-1)<x_u )); % number of upcrossings
    
    y_pk = sqrt(2.0*log(-duration_ratio*Nupcross ./ log(cdf_pk))); % maximum peaks corresponding to specified cumulative probabilities
    CDF_y = stdnormcdf(y_pk);

    % Perform the mapping procedure to compute the CDF of largest peak
    % for X(t) from y(t)
    x_max = stdgaminv(CDF_y,gamma_hi)*beta_hi + mu_hi;
    x_min = stdnorminv(1-CDF_y)*sigma_lo + mu_lo;
    
    % Compute the Mean of the Peaks for process X(t)
    pdf_pk = -y_pk .* cdf_pk .* log(cdf_pk);
    
    if sign(skew_x)>0
        max_est(i) = trapz(y_pk,(pdf_pk .* x_max));
        min_est(i) = trapz(y_pk,(pdf_pk .* x_min));
    else
        max_est(i) = -trapz(y_pk,(pdf_pk .* x_min));
        min_est(i) = -trapz(y_pk,(pdf_pk .* x_max));
    end

end