-
Notifications
You must be signed in to change notification settings - Fork 6
/
Geometric.m
123 lines (101 loc) · 3.74 KB
/
Geometric.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
classdef Geometric < dDiscrete % NWJEFF: Not vectorized
% Geometric(P): Number of trials to first success, including the success (i.e., X=1,2,3,...)
properties(SetAccess = protected)
P, Q
MinPDF % Criterion for deciding when the PDF is small enough to cut off the distribution
end
methods (Static)
function Reals = ParmsToReals(Parms,~)
Reals = NumTrans.Bounded2Real(0,1,Parms(1));
end
function Parms = RealsToParms(Reals,~)
Parms = NumTrans.Real2Bounded(0,1,Reals(1));
end
end
methods
function obj=Geometric(varargin)
obj=obj@dDiscrete('Geometric');
obj.DistType = 'd';
obj.ParmTypes = 'r';
obj.DefaultParmCodes = 'r';
obj.NDistParms = 1;
obj.MinPDF = 1e-10; % was e-20
switch nargin
case 0
case 1
ResetParms(obj,varargin{1});
otherwise
ME = MException('Geometric:Constructor', ...
'Geometric constructor needs 0 or 1 arguments.');
throw(ME);
end
end
function []=ResetParms(obj,newparmvalues)
ClearBeforeResetParmsD(obj);
obj.P = newparmvalues(1);
ReInit(obj);
end
function PerturbParms(obj,ParmCodes)
% Perturb parameter values prior to estimation attempts.
newP = ifelse(ParmCodes(1)=='f', obj.P, obj.P*0.95);
obj.ResetParms(newP);
end
function []=ReInit(obj)
assert(obj.P>0&&obj.P<1,'Geometric P must be 0-1.');
MakeTables(obj);
TrimTables(obj,obj.MinPDF,1-eps(1));
SetBinEdges(obj);
obj.Initialized = true;
if (obj.NameBuilding)
BuildMyName(obj);
end
end
function []=MakeTables(obj)
obj.Q = 1 - obj.P;
obj.LowerBound = 1;
obj.UpperBound = ceil(1+log(obj.MinPDF/obj.P)/log(1-obj.P));
obj.NValues = round(obj.UpperBound - obj.LowerBound) + 1;
obj.DiscreteX = obj.LowerBound:obj.UpperBound;
AllQs = obj.Q*ones(1,obj.NValues);
obj.DiscretePDF = obj.P*cumprod(AllQs)/obj.Q;
obj.DiscreteCDF = cumsum(obj.DiscretePDF);
obj.DiscreteCDF(end) = 1;
end
function thisval=Mean(obj)
if ~obj.Initialized
error(UninitializedError(obj));
end
thisval = 1/obj.P;
end
function thisval=Variance(obj)
if ~obj.Initialized
error(UninitializedError(obj));
end
thisval = (1-obj.P)/obj.P^2;
end
function thisval=Skewness(obj)
if ~obj.Initialized
error(UninitializedError(obj));
end
thisval = (2-obj.P)/sqrt(1-obj.P);
end
function thisval=Kurtosis(obj)
if ~obj.Initialized
error(UninitializedError(obj));
end
thisval = 9 + obj.P^2/(1-obj.P);
end
function s=EstML(obj,Observations,varargin)
if ~obj.Initialized
error(UninitializedError(obj));
end
meanObs = mean(Observations);
ResetParms(obj,1/meanObs);
BuildMyName(obj);
s=obj.StringName;
end
function thisval=Random(obj,varargin)
thisval = geornd(obj.P,varargin{:}) + 1; % geornd is number of failures before first success
end
end % methods
end % class Geometric