The cusp model of catastrophe theory is very closely related to certain multiparameter exponential families of probability density functions. This relationship is exploited to create an estimation theory for the cusp model. An example is presented in which an independent variable has a bifurcation effect on the dependent variable

Cobb, Loren

Munich Personal RePEc Archive

MPRA
Munich Personal RePEc Archive
Estimation Theory for the Cusp
Catastrophe Model
Loren Cobb
University of Colorado Denver
April 1980
Online at https://mpra.ub.uni-muenchen.de/37548/
MPRA Paper No. 37548, posted 22 March 2012 19:16 UTC
Written June 1980 1980 Proceedings of the American Statistical Association 
Published Spring 1981 Section on Survey Research Methods 
Revised June 2010 pp. 772–776 
Estimation Theory for the Cusp Catastrophe Model 
REVISED EDITION1 
Loren Cobb 
Abstract	  
The cusp model of catastrophe theory is very closely related to certain multi-
parameter exponential families of probability density functions. This relationship 
is exploited to create an estimation theory for the cusp model. An example is pre-
sented in which an independent variable has a bifurcation effect on the depend-
ent variable. 
Introduction	  
The elementary catastrophe models of Thom (1975) and Zeeman (1977) have at-
tracted the attention of researchers and theorists throughout the sciences. A per-
sistent problem with virtually all published applications, however, has been the 
absence of statistical procedures for detecting the presence of catastrophes in any 
given body of data. This lack has resulted in some severe criticism of catastrophe 
models being, among other things, speculative and unverifiable (Sussmann and 
Zahler, 1978). Thus catastrophe models have become associated in many minds 
with reckless speculation and intellectual irresponsibility. As part of an effort to 
overcome this problem, this paper presents an estimation theory and the begin-
nings of an inferential theory, in a form useful for survey research applications of 
catastrophe models. 
Catastrophe models come in both dynamic and static forms, the static forms be-
ing simply the equilibria (stable and unstable) of the dynamic forms. The capac-
ity for multiple stable equilibria is inherent in catastrophe models: this is the 
principal feature which distinguishes them from the standard models used in 
linear and polynomial regression. 
In effect, the "control" factors of a catastrophe model correspond to the inde-
pendent variables of a statistical model, and the "behavioral" variable of a catas-
trophe model corresponds to the dependent variable of a statistical model. When 
the control factors are such that the behavioral variable is in a multi-stable situa-
tion, then each stable equilibrium value is a predicted value of the behavioral 
variable — thus there is more than one predicted value. In addition, the unstable 
equilibria which separate the stable equilibria are also predictions of a sort: they 
are the values that we predict that the behavioral variable will not have. This fea-
ture of catastrophe models makes it difficult to define the size of an error of pre-
diction. 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 2 of 11 
There are two ways of overcoming this difficulty. Both of these ways have 
emerged from a study of various forms of dynamic stochastic catastrophe models 
(Cobb, 1978, 1980; Cobb & Watson, 1981). One of these is based on the method of 
moments, while the other is based on maximum likelihood estimation. The latter 
permits hypothesis testing through the use of the chi-square approximation to 
the likelihood ratio test. The former has the advantage of computational simplic-
ity, while the latter is clearly preferable when hypotheses must be tested. 
 
Figure 1: The cusp catastrophe model. 
The	  Cusp	  Model	  
The canonical cusp model can be thought of as a rather peculiar response surface 
model. It's shape may be seen in Figure 1 below. Note that sections taken 
through the depicted surface parallel to the α-axis are just cubic polynomial in y, 
the dependent variable. The entire surface is defined by the implicit equation 
 0 = α + β y − λ
σ
⎛
⎝⎜
⎞
⎠⎟
−
y − λ
σ
⎛
⎝⎜
⎞
⎠⎟
3
. 
If we let z = y − λ
σ
 be the "standardized" dependent variable, then the cubic equa-
tion is simply 
 0 = α + βz − z3 . 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 3 of 11 
It may be seen that λ and σ are location and scale parameters, respectively, for y. 
The roots of the cubic polynomial are the predicted values of z, given α and β. 
When there are three roots, the central root is an "anti-prediction": a prediction of 
where the dependent variable will not be. This feature of the cusp model is clari-
fied in Figure 2, which shows the sequence of conditional probability density 
functions for y, with α fixed as β is increased. This sequence corresponds to the 
trajectory and its projection that are shown in Figure 1. These probability density 
functions will be discussed in a later section. 
 
Figure 2: The cusp family of probability density functions. 
The two dimensions of the control space, α and β, are canonical factors which 
depend upon the actual measured independent variables, say X1,...,Xv . As a first 
approximation, we may suppose that the control factors depend linearly upon 
the independent variables: 
 α = α0 +α1X1 + ...+α vXv ,
β = β0 + β1X1 + ...+ βvXv .
 
Thus the statistical estimation problem is to find estimates for the 2v+4 parame-
ters 
 λ,σ ,α0 ,...,α v ,β0 ,...,βv( ),  
from n observations of the v+1 variables 
 Y ,X1,...,Xv( ).  
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 4 of 11 
As β changes from negative to positive, the conditional probability density func-
tion of y changes in shape from unimodal to bimodal. For this reason the β factor 
will be called the bifurcation factor (it has also been called the splitting factor, by 
Zeeman and others). When α is zero the PDF is symmetrical no matter what the 
value of β. When the PDF is unimodal, α determines its skew: α positive implies 
positive skew, and vice versa. However, when the PDF is bimodal, then α deter-
mines the relative height of the two modes: α positive implies that the right-hand 
mode is higher, and vice versa. To encompass these attributes with a single term, 
α will be called the asymmetry factor (it has also been called the normal factor, a 
rather misleading term in the statistical context). 
Because the model is based on a cubic polynomial, it is possible to define a statis-
tic which discriminates between the unimodal and bimodal cases. This is 
Cardan's discriminant: 
 δ = 12α( )2 − 13 β( )3 .  
When δ is negative the PDF is unimodal, and when it is positive the PDF is bi-
modal. 
Statistical	  Theory	  
The probability density function upon which all of the preceding descriptive sta-
tistics were based is the standard 4-parameter cusp PDF: 
 f (y;α,β,λ,σ ) = ξ exp αz + 12 βz2 − 14 z4( ),   with   z = y − λσ . 
The constant ξ merely normalizes the PDF so that it has unit integral over its 
range, which is the whole real line. The modes and antimodes of the cusp PDF 
may be found by solving ∂f
∂y = 0 . This yields the equation 
 0 = α + βz − z3 , 
which is exactly the same as the implicit equation which defined the cusp surface. 
The modes of the cusp PDF are the predicted values of the cusp model, while the 
antimodes of the cusp PDF are the anti-predictions of the cusp model. The deriva-
tion of the cusp PDF from stochastic catastrophe theory, using stochastic differ-
ential equations, may be found in (Cobb & Watson, 1981). The statistical theory 
was first presented in rudimentary form by (Cobb, 1978). 
The standard cusp PDF can clearly be reparametrized so that it is an exponential 
family, as in: 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 5 of 11 
 f (y) = exp −η + τ1y + τ 2y2 + τ 3y3 + τ 4y4( ).  
Now the well-developed theory (e.g. Lehman, 1959) of exponential families can 
be applied: we know that maximum likelihood estimators (MLEs) exist, are 
unique, and can be found for example by a Newton-Raphson search. This search 
procedure proceeds as follows. Let τ  stand for the vector of parameters 
 τ = τ1,τ 2 ,τ 3,τ 4( ) , 
let S be the vector of sample means defined by 
 Sk = 1n Yik
i=1
n
∑ ,   for  k = 1,2, 3,4,  
let Mk (τ )  be the vector of expectations of Y given τ defined by 
 Mk (τ ) = E Y k⎡⎣ ⎤⎦,   for  k = 1,2,3,4,  
and let H (τ ) be the 4x4 covariance matrix defined by 
 Hij (τ ) = Cov Y i ,Y j⎡⎣ ⎤⎦,   for  i, j = 1,2,3,4.  
The Newton-Raphson search is an iterative procedure which starts from an ini-
tial guess for τ, say t0 . Each iteration is given by the vector expression 
 tk+1 = tk + S − M (tk )( )H (tk )−1.  
This calculation is performed repetitively until S = M (t) , within the limits of 
computational accuracy. 
It should be noted that after each iteration the vector M(t) and the matrix H(t) 
must be recalculated. These moments must be found by numerical integration, 
since closed-form expressions for the moments of the cusp PDF are not known.2 
The preceding discussion applies to the estimation of the four parameters of the 
cusp PDF, given observations of the variable Y. If, however, the cusp PDF is to be 
used as the conditional density of Y, given the values of the independent vari-
ables, then the maximum likelihood procedure becomes more complicated. If we 
use the previously-stated assumption that the factors α and β are linear combina-
tions of the independent variables, then the extension of the Newton-Raphson 
technique for finding the maximum likelihood estimates is straightforward. 
There are now 2v+4 parameters to be estimated, and the matrix H(t) becomes 
(2v + 4) × (2v + 4)  dimensional. 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 6 of 11 
We now proceed to estimation by the method of moments. It will be seen that 
this method, in contrast to the method of maximum likelihood, is extremely easy 
to implement. 
Even though closed-form expressions for the moments of the cusp PDF do not 
exist, moment estimators are trivial to derive with the aid of the following gen-
eral theorem: 
Theorem:  Let g(x, y)  be a polynomial function of x and y such that 
 0 < exp − g(x, y)dy∫⎡⎣ ⎤⎦
−∞
+∞
∫ < ∞,   ∀x.  
Let ψ (x)  be the reciprocal of this quantity. Suppose that a random variable Y de-
pends upon x in such a way that its conditional density is given by 
 f (y x) =ψ (x)exp − g(x, y)dy∫⎡⎣ ⎤⎦.  
Assume that the joint density of X and Y has moments of all orders, and let h(x)  
denote the density of the random variable X. Then for any non-negative j and k, 
 E X jY kg X,Y( )⎡⎣ ⎤⎦ = kE X jY k−1⎡⎣ ⎤⎦.  
Proof:  Note that f (y x)  is asymptotically zero as y tends to either +∞  or −∞ . 
Since g(x, y)  is a polynomial, we also have that yk f (y x)  to zero in the same way. 
Further, we can write g(x, y)  as 
 g(x, y) = − ∂
∂y ln f (y x( ) = −
∂f (y x) ∂y
f (y x) . 
Substituting this expression into the moment formula, we obtain 
 
E X jY kg X,Y( )⎡⎣ ⎤⎦ = x jykg(x, y) f (y x)h(x)dydx∫∫
= x jyk −∂f (y x)( )h(x)dx∫∫
= x jh(x) − yk∂f (y x)∫( )dx∫ .
 
Now use integration by parts on the inner integral, and observe that one of the 
parts is identically zero: 
 
− yk∂f (y x)∫ = − yk f (y x) −∞
+∞
+ k yk−1 f (y x)dy∫
= 0 + k yk−1 f (y x)dy∫ .
 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 7 of 11 
Thus we now have 
 
E X jY kg X,Y( )⎡⎣ ⎤⎦ = x jh(x) k yk−1 f (y | x)dy∫( )dx∫
= k x jyk−1h(x)dydx∫∫
= kE X jY k−1⎡⎣ ⎤⎦,
 
which was to be shown.  
This theorem enables the method of moments to be applied to models that, like 
the elementary catastrophes, are expressed as implicit equations. Before examin-
ing the cusp model, it may be worthwhile to show how it can be applied to ordi-
nary linear regression. The linear regression model 
 y = a + bx + ε  
can be written in implicit equation form as 
 g(x, y) = y − a − bx
σ 2
= 0,  
where σ 2  will turn out to be the variance of ε . The conditional PDF of y given x 
is 
 
f (y x) =ψ (x)exp −
1
2 y2 − ay − bxy
σ 2
⎡
⎣
⎢
⎤
⎦
⎥
=ψ 1(x)exp − 12
y − (a + bx)
σ
⎛
⎝⎜
⎞
⎠⎟
2⎡
⎣
⎢
⎢
⎤
⎦
⎥
⎥
.
 
This is clearly a parametrized family of normal densities, N a + bx,σ 2( ) . (To ob-
tain this formula, simply complete the square and absorb the terms in x into the 
function ψ, the normalizing constant.) 
To find estimation equations for a and b, use the theorem twice, first with j = k = 0 
and second with j = 1 and k = 0: 
1. E g X,Y( )⎡⎣ ⎤⎦ = 0    ⇒     a + bE X[ ] = E Y[ ],  
2. E Xg X,Y( )⎡⎣ ⎤⎦ = 0    ⇒     aE X[ ] + bE X 2⎡⎣ ⎤⎦ = E XY[ ],  
Notice that when sample moments are substituted for these expectations, we ob-
tain the usual Gauss-Markov normal equations for linear regression. 
To estimate σ 2 , use the theorem again, this time with j = 0 and k = 1: 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 8 of 11 
3. E Yg X,Y( )⎡⎣ ⎤⎦ = 1    ⇒     E Y 2⎡⎣ ⎤⎦ − aE Y[ ] + bE XY[ ] = σ 2 ,  
which is the correct formula for the residual variance of Y after the linear effect of 
X has been removed by linear regression. 
Turning now to the cusp model, let us consider first the model with no inde-
pendent variables: 
 g(y) = a + by + cy2 + dy3,     d > 0( ).  
This model has a PDF given by 
 f (y) = ξ exp ay + 12 by2 + 13 cy3 + 14 dy4⎡⎣ ⎤⎦,  
which has modes and antimodes at the roots of g(y) = 0 . The transformation of 
the parameters from a,b,c,d( )  to the standard coefficients α,β,λ,σ( )  is accom-
plished by: 
 
σ = d4 ,
λ = − c3d ,
α = −σ a + bλ + cλ2 + dλ 3( ),    and
β = −σ 2 b + cλ( ).
 
Estimation of the parameter vector a,b,c,d( )  proceeds from an application of the 
theorem. Let 
 µk = E Y k⎡⎣ ⎤⎦.  
From a single application of the theorem we can derive a linear difference equa-
tion (with on varying coefficient) for the moments of Y: 
 E Y
kg(Y )⎡⎣ ⎤⎦ = kE Y k−1⎡⎣ ⎤⎦
⇒ kµk−1 = aµk + bµk+1 + cµk+2 + dµk+3.
 
Simply apply this result with  k = 0, 1, 2, 3  to obtain a system of four linear equa-
tions in the four unknowns a,b,c,d( ) . Substitute sample moments for the expec-
tations, and solve the system. Transform the resulting estimates as indicated 
above to obtain α,β,λ,σ( ) . 
It is trivial to extend this technique to models with independent variables. For 
example, suppose there is one independent variable, say X. Then the model is 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 9 of 11 
 g(x, y) = b1 + b2x + b3y + b4xy + b5y2 + b6y3. 
Estimation of the coefficients b1,b2 ,b3,b4 ,b5 ,b6( )  proceeds as before. Apply the 
theorem six times, with j = 0 and k = 0, 1, 2, 3, and then with j = 1 and k = 0, 1. 
Substitute sample moments for expectations, and solve the resulting system. Fi-
nally, the standard coefficients are obtained from 
 
σ = b64 ,
λ = − b53b6
,
α0 = −σ b1 + b3λ + b5λ2 + b6λ 3( ),
α1 = −σ b2 + b4λ( ),
β0 = −σ
2 b3 + b5λ( ),
β1 = −σ
2b4 .
 
Estimation by the method of moments does not yield estimators with known 
sampling distributions, and cannot easily be used for hypothesis testing [Note: 
The approximate sampling distribution of these estimators was later derived by Ferdon in 
her PhD dissertation (Ferdon, 1983).] The maximum likelihood method yields esti-
mators that are efficient and that have known asymptotic sampling distributions, 
but they are, from a computational point of view, inefficient. It is possible, how-
ever, to use the moment estimates as the initial guess for the Newton-Raphson 
iterations, thus cutting down somewhat the time required to calculate the MLEs. 
Once MLEs have been obtained, it is possible to test hypotheses — for example, 
to compare a linear regression model to the equivalent cusp catastrophe model 
— using the chi-square approximation to the likelihood ratio test. 
An	  Example	  
An excellent example of published empirical data which seems to exhibit a bifur-
cation in the dependent variable has been quoted by (Zeeman, 1977, pp. 373–85). 
The data come from a study of driving performance before and after the inges-
tion of alcohol (Drew, Colquhoun, and Long, 1959). Essentially, the authors 
found that the change in time per lap (i.e. the driving speed) was strongly af-
fected by the position of the subjects on the Bernreuter scale of introversion. 
However, as is visible from Figure 3, it is clear that whereas extroverts continued 
to drive at about the same speed after drinking, the introverts either drove faster 
or slower, and few stayed at the same speed. 
These data were reproduced as a figure in (Zeeman, 1977, Fig. 1), from which 
approximate data were recovered by digitization. Following Zeeman, three cases 
were eliminated as extreme outliers, leaving the 37 cases depicted in Figure 3. 
The six-parameter cusp model with one independent variable was fitted to the 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 10 of 11 
data using the method of moments as given above, and the resulting relationship 
between change in driving speed after alcohol (the dependent variable) and in-
troversion (the independent variable) is shown. The dashed line indicates values 
that are predicted not to occur (the anti-predictions). 
 
Figure 3: A bifurcation model fitted to data. 
Zeeman also used a cusp model in his article, although it differs substantially 
from the one in Figure 3. Poston and Stewart (1978, pp. 420–23) criticized Zee-
man's model on psychological grounds, and suggested the bifurcation model 
that appears here. 
References	  
Cobb, Loren (1978) “Stochastic catastrophe models and multimodal distribu-
tions,” Behavioral Science, v. 23, 360–374. 
Cobb, Loren (1980) “Stochastic differential equations for the social sciences,” 
Chapter 2 in Mathematical Frontiers of the Social and Policy Sciences, edited 
by Loren Cobb and R. M. Thrall. Boulder, CO: Westview Press. 
Cobb, Loren and Watson, William B. (1981) “Statistical catastrophe theory: An 
overview,” Mathematical Modeling, v. 1, 311–317. 
Drew, G. C., Colquhoun, W.P., and Long, H.A. (1959) “Effect of small doses of 
alcohol on a skill resembling driving,” Memo 38. London: Medical Re-
search Council. 
Loren Cobb Estimation Theory for the Cusp Catastrophe Model Page 11 of 11 
Ferdon, Marybeth E. (1983) Inference for Quadratic and Catastrophe Response Surface 
Models. Unpublished doctoral dissertation. Charleston, SC: Department of 
Biostatistics and Epidemiology, Medical University of South Carolina. 
Grasman, Raoul P., van der Maas, Han L., and Wagenmakers, Eric-Jan (2009) 
“Fitting the cusp catastrophe in R: A cusp package primer,” Journal of Sta-
tistical Software, vol. 32, #8, 1–27. 
Lehman, E.L. (1959) Testing Statistical Hypotheses. New York: Wiley. 
Poston, Timothy and Stewart, Ian (1978) Catastrophe Theory and Its Applications. 
London: Pitman. 
Sussmann, H.J. and Zahler, R.S. (1978) “A critique of applied catastrophe theory 
in the behavioral sciences,” Behavioral Science, vol. 23, 383–389. 
Thom, René (1975) Structural Stability and Morphogenesis. Reading, MA: Benjamin. 
Zeeman, E. Christopher (1977) Catastrophe Theory: Selected Papers. Reading, MA: 
Addison-Wesley. 
Endnotes	  
                                                
1 This non-refereed proceedings paper is one of my most-cited articles. For this revision I have 
fixed several grammatical and mathematical errors, and brought some references up to date. 
2 There are efficient alternatives to numerical integration, as employed for example in the R pack-
age for cusp estimation by (Grasman et al, 2009). 


Estimation Theory for the Cusp Catastrophe Model

https://mpra.ub.uni-muenchen.de/37548/1/MPRA_paper_37548.pdf

Estimation Theory for the Cusp Catastrophe Model

Abstract

Similar works

Full text

Available Versions

Munich Personal RePEc Archive