WGGAMSTAT Mean and variance for the Generalized Gamma distribution. CALL: [m,v] = wggamstat(a,b,c) m, v = the mean and variance, respectively a,b,c = parameters of the Generalized Gamma distribution. (see wggampdf) Mean (m) and variance (v) for the Generalized Gamma distribution is m=c*gamma(a+1/b)/gamma(a) and v=c^2*(gamma(a+2/b)/gamma(a)-gamma(a+1/b)^2/gamma(a)^2); Example: param = {1,2,4};N = 1000; x = wggamrnd(param{:},N,1); [mean(x),var(x)] [m,v] = wggamstat(param{:})
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Gamma function. | |
Not-a-Number. |
001 function [m,v]= wggamstat(a,b,c); 002 %WGGAMSTAT Mean and variance for the Generalized Gamma distribution. 003 % 004 % CALL: [m,v] = wggamstat(a,b,c) 005 % 006 % m, v = the mean and variance, respectively 007 % a,b,c = parameters of the Generalized Gamma distribution. (see wggampdf) 008 % 009 % Mean (m) and variance (v) for the Generalized Gamma distribution is 010 % 011 % m=c*gamma(a+1/b)/gamma(a) and 012 % v=c^2*(gamma(a+2/b)/gamma(a)-gamma(a+1/b)^2/gamma(a)^2); 013 % 014 % Example: 015 % param = {1,2,4};N = 1000; 016 % x = wggamrnd(param{:},N,1); 017 % [mean(x),var(x)] 018 % [m,v] = wggamstat(param{:}) 019 020 021 % Reference: Cohen & Whittle, (1988) "Parameter Estimation in Reliability 022 % and Life Span Models", p. 220 ff, Marcel Dekker. 023 024 % Tested on; Matlab 5.3 025 % History: 026 % revised pab Jul2004 027 % fixed a bug 028 % reviseed pab Dec2003 029 % fixed abug: k1 -> k3 030 % revised pab 24.10.2000 031 % - added comnsize, nargchk + default value for b and c 032 % added ms 09.08.2000 033 error(nargchk(1,3,nargin)) 034 if nargin<2, b=1;end 035 if nargin<3, c=1;end 036 [errorcode a b c] = comnsize(a,b,c); 037 if errorcode > 0 038 error('a b and c must be of common size or scalar.'); 039 end 040 041 % Initialize m and v to zero. 042 m = zeros(size(a)); 043 v=m; 044 ok = (a > 0 & b>0 & c>0); 045 k=find(ok); 046 if any(k), 047 m(k) = b(k).*gamma(a(k)+1./c(k))./gamma(a(k)); 048 v(k) = b(k).^2.*gamma(a(k)+2./c(k))./gamma(a(k))-m(k).^2; 049 end 050 051 k3 = find(~ok); 052 if any(k3) 053 tmp = NaN; 054 m(k3) = tmp(ones(size(k3))); 055 v(k3)=m(k3); 056 end 057 058 059 060 061 062
Comments or corrections to the WAFO group