Experiments in Nonlinear Psychophysics:
Replications of Some of Robert Gregson's Work

Psychophysics, the study of human responses to physical stimulation, was dominated first by Fechner's logarithmic law, and then by Stevens' power law. Both laws assume that changes in responses within the organism are directly proportional to changes in the level of external stimulation. While useful in certain contexts, these laws do not do justice to the full complexity of human psychophysical response; and many researchers have pointed out the need for a more sophisticated mathematical psychophysics.

One promising direction, in this regard, is nonlinear science. In recent studies of the brain, neuroscientists have begun to realize that perceptual processes often reflect patterns found in non-linear dynamics. Walter Freeman, for example, has shown that some apparently random perceptual processes actually have deterministic dynamics, with chaotic attractors. In 1991, he published a paper "The physiology of perception," reporting experiments in which he trained rabbits to recognize several different smells. He attached electrodes to the olfactory bulb surface of the rabbits in order to collect sets of 60 to 64 simultaneously recorded electroencephalogram (EEG) tracings. In the experiments, the EEG tracings from the electrodes are "as unpredictable and irregular as freehand scrawls" [1]. Yet they manifest perceptual information. This can be shown by "phase portraits" which are made from the EEGs generated by a computer model of the brain. In Freeman's experiments, the phase portraits reflect the overall activity of the olfactory at rest and during perception of a familiar scent. As he states, "Resemblance of the portraits to irregularly shaped, but still structured, coils of wire reveals that brain activity in both conditions is chaotic: complex but having some underlying order." He also indicates that the olfactory EEGs are closer to being periodic during perception than during rest.

Recognition of the complexity and nonlinearity of perception has led some researchers to propose the use of non-linear dynamics models for describing the perceptual processes. Most impressively, Robert Gregson has developed a non-linear psychophysical function he calls the "gamma recursion" and has conducted many computer simulations of psychophysical phenomena using the model. The model has been shown to demonstrate many complex, psychophysically plausible behaviors [2].

Gregson's argument for the non-linearity of the perceptual processes is as follows: "It is known that sensory pathways at higher cortical levels can be multiple, recursive, perseverative, and nonlinear. There are also grounds for treating large aggregates of connected neural structures as capable of generating nonlinear and chaotic dynamics as part of their normal functioning. Given those constraints, a cubic complex polynomial recursion gamma was defined, which satisfies the criteria stipulated. This model generates, intrinsically, both observable response and unobservable (theoretical, deterministic) noise at the same time, without excluding the fact that additional stochastic noise can also exist in the system modeled" [3].

This paper reports an attempt to replicate several of Gregson's results. In many cases my specific numerical results do not agree with Gregson's; however, the basic qualitative pattern of the results is the same. In accordance with Gregson, I conclude that the gamma recursion is a viable tool for studying the nonlinearity underlying human psychophysical responses.

The "original" (one-dimensional) gamma recursion is introduced for the first time in 1988 in Gregson's book Nonlinear Psychophysical Dynamics. It is defined as

Y_j+1 = -a (Y_j - 1)(Y_j + ie)(Y_j - ie), i² = -1

where a is real (Re), ie is imaginary (Im), Y is complex (Re, Im), and the boundary conditions are 0 < Y(Re) < 1, Y₀ = (), .5, < 10^-8 (i.e. the initial condition .5 + 10^-9 i), 2 < a < 4, 0 < e < .5, and very approximately ae < 1.7 in the region of .5 < e < .7 to avoid the explosive condition [2].

An interesting point of this function is that with a little modification of the parameter value a, the gamma recursion can model the psychometric function that is a CNO-like function. For example, one of the modification in parameter value a of the gamma, he calls [[Gamma]]V7, is defined as

0 < U < 1, a_min < a < 4, e fixed

(number of iterations) fixed, Y₀ = .5 , [[epsilon]] (complex)

a = a_min + (3.99 - a_min) U,

where U is the stimulus, and Y(Re) is the response [2]. The shapes of these functions after 10 iterations ( = 10) varies completely as e changes from .03 to .3. Figure 1 shows the change over that range. And with a_min = 3.2 and e = .15, the function shows the pattern of a psychometric function.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 Curves produced by V7 and a CNO curve

Figure 1 (a)(e) shows the curves produced with the real part of the gamma recursion after 10 iterations with a = 3.2 + .79 U (a_min = 3.2), 0 < U < 1. (a) e = .03,

(b) e = .1, (c) e = .15 (shows a CNO-like curve), (d) e = .2, (e) e = .3.

Figure 1 (f) shows the cumulative normal ogive, which is a psychometric function generated by psychophysical data.

The x-axis represents the stimuli, and the y-axis represent the responses.

The values on the graphs are scaled for graphing purposes.

The gamma recursion models the transformation from stimulus to response in a single channel. However, to represent, for example, bilateral stimuli inputs such as fully vision or audition requires two dimensional inputs. For this purpose, Gregson introduced the two-channel gamma model he calls the 2[[Gamma]]. The system 2[[Gamma]] is classified into two cases by the difference of their cross-coupling between channels; the first case is simpler, however, it would not generate the irregularity which might, for example, characterize some visual illusions [5]. The 2[[Gamma]] can be extended to more general vectorial and lattice models called n and (nn).

The descriptions of Case 1 and Case 2, and the model are extracted from Gregson's book n-Dimensional Nonlinear Psychophysics:

Case 1: If the (Im) output of Y_obs, (which is the value of Y after recursions), is taken as interpretable as a form of internal noise, as has been shown to be plausible in other contexts (for example, it exhibits intermittency), then take the total within-loop noise at any step in the recursion to be the union of the component noises. So in each cycle the two Y(Im) are replaced by the term max(Y_1,j, Y_2,j)(Im) and as the (Im) components go into limit cycling more readily (that is, at lower system parameter values) than their (Re) counterparts, the effect will depend crucially on the set of values (a₁, e₁, a₂, e₂) (Campbell and Gregson, 1990). It can be shown, by simulation, that high (e₁, e₂) values induce some curious effects for large U_j (inputs); a paradoxical "turn-over" effect in

U Y_obs graphs for at least one component is observable (Gregson and Britton, 1990). At the same time the other component may create a monotonically decreasing graph of U Y_obs.

For a double loop put, for case 1, where the two dimensions are denoted by h and i suffices,

Y_hij(Im) = max(Y_hj(Im), Y_ij(Im))

and then define

Y^*_hj = (Y_hj(Re), Y_hij(Im))

and

Y^*_ij = (Y_ij(Re), Y_hij(Im).

Case 2: A second cross-linkage is possible if the e terms are made to be functions of the a terms in the opposite dimension. This represents a process of mutual inhibition. Each dimension is desensitized to differences in its own input as the input to the other dimension is increased, viz.:

e_h = _h a_i^-1, e_i = _i a_h^-1,

where _h > 0, _i > 0, and 2.6 < a_min < 4.4 is a working range for 0.05 < e < 0.40, and approximately .20 < < 1.2 will avoid the explosive condition of ae > 1.7.

The advantage of writing e in terms of and a is that the number of parameters in the model can be reduced by one; e₁ and e₂ are removed and is introduced, with

_h = _i = . For case 2, additionally we have

Y_h(j+1) = -a_h (Y^*_hj - 1)( Y^*_hj - i_h (a_i)^-1)( Y^*_hj + i_h (a_i)^-1)

and correspondingly

Y_i(j+1) = -a_i (Y^*_ij - 1)( Y^*_ij - i_i (a_h)^-1)( Y^*_ij + i_i (a_h)^-1).

See figure 2 for better understanding of the model.

Pure Delay Filter

Sensory ₁

U₁ Dimension a₁ e₁ Y_obs,1

One

Inputs Real Observable

Sensory ₂

U₂ Dimension a₂ e₂ Y_obs,2

Two

Figure 2 Diagram of cross-coupling used in 2[[Gamma]] case2, the labelled ₁, ₂ are themselves the recursions.

To show non-linear psychophysical dynamics using his 2[[Gamma]] , Gregson has carried out several studies with psychophysical data. In one study, Gregson [3] applied 2[[Gamma]] to approach Fechner's Paradox. Fechner's Paradox is where apparent brightness of an object viewed binocularly can, under conditions where the input to one eye is diminished by filtering, be less than its brightness viewed monocularly by the unfiltered eye. The same phenomenon in pooling two sensory inputs has its analogues in audition and in olfaction. In this case also, Gregson points out that the analysis with a linear model makes a limited success; "The idea that unidimensional stimuli coalesce vectorialy to generate mixtures of two stimuli, in terms of intensities, has been used with limited success in modelling various sensory modalities; in the context of Fechner's Paradox it was advanced by Curtis and Rule (1978). It is now known that it fails in vision (Irtel, 1986), as it also fails in olfaction (Gregson, 1986)" [3]. Using his 2[[Gamma]], case2, however, Gregson states that "great variations in the form and scatter of isointensity (brightness or loudness) distributions occur in the region of a₁/a₂ [two input stimuli] and a model should predict such local instability, as well as predicting much clearer trends in the midrange where a₁ and a₂ approach equally" [3]. Gregson's results indeed predict the instability.

I replicated Gregson's research using the "original" gamma recursion with the same parameter values of a and e. While I obtained several results that were consistent with Gregson's, there were a number of significant discrepancies. Figure 3 shows the chaotic behavior of the gamma recursion with certain parameter values a, and e, where my results support Gregson's findings. Figure 4 shows result of the gamma recursion for some parameter values a and e where Gregson's results are not supported.

(a)

(b)

Figure 3 Phase portraits showing chaotic behavior in the real and imaginary domain for the gamma recursion with specific parameter values of a and e.

(a) parameter values a = 5.75 and e = 0.2

(b) parameter values a = 4.20 and e = 0.4340

Each point is plotted as an x-y coordinate where x-axis is the real/imaginary value of Y at the jth iteration and y-axis is the real/imaginary value of Y at the j+1th iteration except the third one of (b). In this case, x-axis is the real value of Y at the jth iteration and y-axis is the imaginary value of Y also at the jth iteration.) .

Each successive point is joined, providing a trace of the output of the gamma recursion across iterations. Complex traces reflect chaotic behavior. The values after each iteration are unique and seem random. However, they are confined within a certain region forming a pattern that is called a "strange attractor".

Figure 4 Phase portraits showing the dynamics in the real and imaginary domain for the gamma recursion with specific parameter values of a and e.

Comparison is made between Gregson's results (left) and mine (right).

(a) parameter values a = 3.31 and e = 0.2. With these parameter values, the gamma recursion converges to 0.337403 + 4.59355 10^-16i. Therefore Gregson's result is not accurate.

(b) parameter values a = 6.543 and e =0.0. The result I obtained (right) shows the trace of up to 29 iterations. At the 30th iteration, the real value of Y escapes the region of (-1,1), and therefore causes the explosion.

In order to probe the characteristics of the one dimensional gamma recursion further, Gregson also conducted an experiment to see how the Y_obs (observation, the real component of Y) varies over the parameter values a and e for fixed recursions. His reasoning for this action is as follows;

"If [2.2] [the gamma recursion] runs a limit of recursions then the result is Y_obs, which as defined is complex. It may be single-valued or oscillating in one or both components. For a fixed , over the parameter space {a, e} we may thus construct a bivariate response surface P_Y, of { Y_obs, a, e}, and this at any point will have four values; minY(Re), maxY(Re), minY(Im), maxY(Im). There are thus in theory four distinct response surfaces, but they are not distinguishable at all (a, e) points, and not all of equal interest" [2].

Figure 5 shows Gregson's results of P_Y, at = 1000 and = 1001 which are the surface graphs of Y at 1000th iteration and 1001st iteration over the region of

2.5 < a < 6.5 and 0.0 < e < 0.8. Figure 6 shows my results with same conditions as Gregson's. My results are clearly different to those of Gregson's.

---

(a)

(b)

Figure 5 Gregson's results for the gamma recursion after 1000 (a) and 1001 (b) iterations with parameter values over the region of 2.5 < a < 6.5 and 0.0 < e < 0.8.

(a)

(b)

Figure 6 My results for the gamma recursion after 1000 (a) and 1001 (b) iterations with parameter values over the region of 2.5 < a < 6.5 and 0.0 < e < 0.8.

In addition to replicating experiments reported in Gregson's books, I also reviewed Gregson's unpublished paper "The second eigenvalue in 2 gamma case two dynamics and dimensionality of output". In the paper, Gregson examines the dynamics of a local region in 2[[Gamma]] by utilizing the instability of the second (smaller) eigenvalue of the system. Then he applies the theory to Eisler cusps, a peculiar phenomenon observed in the task of estimating dimensions of rectangles. As Gregson observes, "When subjective width with varying heights are plotted against shape (ratio of height over width) with physical width as parameter, curves with cusps are obtained" (Eisler, Eisler and Gregson, 1995). It has been named a `cusp'" [6].

In his study of the Eisler cusp, Gregson studies the eigenvalues of 2[[Gamma]] Case 2 by constructing a 2 x 2 matrix of the form

Y_2,j =