A novel grey–fuzzy–Markov and pattern recognition model for industrial accident forecasting

Industrial forecasting is a top-echelon research domain, which has over the past several years experienced highly provocative research discussions. The scope of this research domain continues to expand due to the continuous knowledge ignition motivated by scholars in the area. So, more intelligent and intellectual contributions on current research issues in the accident domain will potentially spark more lively academic, value-added discussions that will be of practical significance to members of the safety community. In this communication, a new grey–fuzzy–Markov time series model, developed from nondifferential grey interval analytical framework has been presented for the first time. This instrument forecasts future accident occurrences under timeinvariance assumption. The actual contribution made in the article is to recognise accident occurrence patterns and decompose them into grey state principal pattern components. The architectural framework of the developed grey–fuzzy– Markov pattern recognition (GFMAPR) model has four stages: fuzzification, smoothening, defuzzification and whitenisation. The results of application of the developed novel model signify that forecasting could be effectively carried out under uncertain conditions and hence, positions the model as a distinctly superior tool for accident forecasting investigations. The novelty of thework lies in the capability of themodel inmaking highly accurate predictions and forecasts based on the availability of small or incomplete accident data.


Introduction
An industrial accident refers to an undesirable, unanticipated and uncontrollable event potentially capable of producing injuries, losses of lives, asset destruction, and disturbance to social as well as economic activities or even leading to degradation of the environment in an industrial system. An accident is an occurrence triggered by human or non-human (i.e. entities, materials or emissions) in which the worker engaged in service to the industry may be injured. Every year, numerous literature reports are given, which declare an increasing number of industrial accidents globally. As a result of concerns to control accident occurrences, accident investigations are now a vital part of scientific reporting and a requirement by government agencies to all industrial organisations worldwide. Government policy stipulates proper reporting of accidents, its control and management. Hence globally, industrial managers are taking advantage of sound scientific studies to adopt models for their industries bearing in mind that an improperly planned accident control scheme could lead to substantial monetary losses due to accident claims.
For some years now, employing industrial forecasting models in accident forecasting relying on multiple factors has been justified by the fact that causal factors of accidents are attributed to human, equipment and managerial deficiencies (Cooke and Rohleder, 2006;Mohaghegh et al. 2009;Rathnayaka et al. 2011). Although proponents of models further justify the use of multiple causal factors, they also acknowledge that their degrees of interactions are also complex (Qureshi 2008;Stringfellow 2010). Unfortunately, since the institution of multivariate scales in the prediction and forecasting of industrial accidents, there has been a broad-spectrum of criticisms regarding the existence of unsatisfactory results. The problems of multivariate models are due to (1) inability to thoroughly capture the levels of interactions; (2) uncertainties; (3) randomness (Mao and Sun 2011);and (4) imprecision (Zheng and Liu 2009) inherent in accident causes and occurrences. In the sense of solving the problems attributed to multivariate models, the classical univariate prediction models (UPMs) were developed. UPMs are classical predictive models such as auto-regression and integrated moving average (ARIMA), exponential smoothing (ESM) and moving average (MA) adapted and applied in industrial accident forecasting (Kim et al. 2011;Kang et al. 2012;Aidoo and Eshun 2012). Scholars, however, also significantly criticised UPMs, in recent times. According to the literature, these mentioned models have not been entirely accurate in their applications to forecasting industrial accident occurrences. The drawbacks attributed to these models are as follows: take the case of MA, a constant mean of occurrence is assumed. However, this may not be true in practical instances of real-time occurrences. Another weakness of UPMs may be picked from the ARIMA model. It requires the availability of extensive data sizes to be able to make dependable predictions (Brockwell and Davis 2002). But this requirement of a large data size in the industrial world characterised by rapid information changes is a luxury that may be difficult to attain. In addition, accessing information in less industrially developed economies is quite challenging and such models may not be applicable in such environments.
Thus, for the aforementioned issues, we consider both UPMs and multivariate prediction models inappropriate for industrial accident forecasting. Yet there must be progress in the field. As the world experiences breakthrough in research on soft computing tool, more areas in science and technology are adopting these tools in their areas. Therefore, more recently, there has been a huge shift in focus towards accident occurrence prediction using non-traditional artificial intelligence (NTAI) forecasting approaches. With NTAI, new knowledge frontiers have been given birth to, expected to radically explode to benefit members of the industrial accident community. Models such as the artificial neural network (ANN) (Zheng and Liu 2009;Oraee et al. 2011), genetic algorithm (GA) (Farahat and Talaat 2012); grey (GM) (Jiang 2007;Lan and Ying 2014), grey-Markov model (Zhang 2010;Mao and Sun 2011;Huang et al. 2012a, b) and fuzzy time series models (FTSMS) (Khev and Yerpude 2015) have been employed in their original or modified forms for forecasting accidents which occur during mining, construction, transportation and processing activities. The results obtained from these model applications in industrial accident forecasting have also been very encouraging.
The organisation of the current work is as follows: the motivation and study objectives are stated in Sect. 1. A review of the literature is presented in Sect. 2. Sections 3, 4, 5, 6, and 7 are devoted to discussing the methodology of the proposed model. Model tests and validation results and discussion are given in Sect. 8. Conclusions concerning the model are shown in Sect. 9 alongside related future research directions.

Related literature
The application of grey, fuzzy and Markov principles in forecasting, as single concepts or merged together in different combination formats has begun to gain increasing popularity in recent years. In this section, a review of literature is given. Grey-fuzzy-Markov (GFM) forecasting technique is a hybrid model which combines the characteristics of the grey, fuzzy and Markov models. GFM models have been developed based on the understanding that hybrid models have greater forecasting potentials than single evaluation models (Li and Li 2015). GFMs have found applications in areas such as electrical load analysis (Asrari et al. 2012) and biofuel production (Geng et al. 2015). The grey aspect of the GFM has its major focus on uncertainty inherent in sparsely available information (Deng 1982;Liu 2011). The model has been deeply explored for forecasting purposes and is evident by the development of several forms of it. Generally, a grey system can be mathematically expressed as a 0AE 2 a _ ; â h i ð1Þ a 0AE is a crisp value or an interval and exist as a component of a base set or interval ½a _ ; â. Basic arithmetic, properties such as addition and multiplication as well as associative and commutative properties also apply in grey systems analysis (Hickey et al. 2001;Arroyo et al. 2011).
Two general forms of grey models, namely differential transfer function-based models (DTFM) and interval arithmetic-based models (IAM) have been mainly employed in areas such as energy consumption, finance and equipment degradation for crisp value forecasting (Kayacan et al. 2010;Tangkuman and Yang 2011;Mostafaei and Kardooni 2012) and interval forecasting (Garcia-Ascanio and Mate 2010; Zhao et al. 2014), respectively. DTFM involves the use of sequence operators ) and unique mathematical representations to describe inputs and outputs under the assumption of exponential data behaviour. The most popularly employed grey model is the GM(1,1). It has been used in its pure form (Zhao et al. 2014;Tong 2016) or modified forms (Mao and Chirwa 2006;Jiang 2007;Zhang 2010;Mao and Sun 2011) for accident forecasting. The IAM involves the application of mathematical operations on grey intervals created from data to produce degeneration or interval forecast. The presence of the grey component, GFM, enables it to make accurate forecasts in the presence of limited and incomplete data, the fuzzy component of the model functions to eliminate the problem of vagueness and uncertainty in data (Chen and Hsu 2004;Kher and Yerpude 2015), while the Markov component deals with problems concerning fluctuating and random occurrences (Geng et al. 2015).
Generally, the procedure of GFM forecasting using the GM(1,1) as part of its component involves three major stages (Huang et al. 2012a, b;Li and Li 2015).
Stage 1: Building the grey model This involves: 1. The creation of a time sequence for the set of available collection of industrial accident data 2. Passage of the created sequence into an accumulated generating operation (AGO). A modified sequence is obtained in the process 4. Solving to obtain the grey parameters a and b, GM(1,1) forecast is then obtained as Stage 2: Fuzzy classification of grey model errors. This involves the linguistic classification of the percentage errors e k of each model forecast into j number of classes carried out under the assumption of time invariance data behaviour (Sullivan and Woodward 1994). By the use of membership functions, the membership of e k in each fuzzy class m lðe k ; mÞ; m : 1:2:3; . . .; j ½ is established. Huang et al. (2012a, b) and Li and Li (2015) employed the maximum membership principle max lðe k ; mÞ ½ to establish the actual fuzzy class in which e k belongs.
Stage 3: Markov state transition On the assumption that m m : 1; 2; 3; . . .; j ð Þexist as a Markov chain of states s m bounded by ðs mL ; s mU Þ, a Markov transition matrix which shows the probability of transition of the state in which e k belongs is sðe k Þ, from its current state y to another state z ðP yz Þ in t stepwise time order is set up: 3. The GFM forecast for time step t is finally obtained as Using this technique, Huang et al. (2012a, b) employed a dynamic grey model in detecting the dynamic trend of accident fatalities in the construction industry. Li and Li (2015) used an unbiased GM (1,1) based GMF in also forecasting construction accidents.
This technique has been shown to improve on GM(1,1) and grey-Markov model prediction accuracies. However, the degree of prediction accuracies is limited. This is because the technique is actually directed at grey model prediction correction and as such, their prediction accuracies are directly dependent on the prediction accuracy of the grey model base. Thus, situations may exist in which GFM variants may not make be able to make forecasts that show significant improvement over those of the GM(1,1) base component.
In addition, Markov-chain transition analysis using the classical Markov state probability matrices and relations only provide general information on data dynamics. This is because the approach requires the availability of specific pre-existing states having similar characteristics to the current state occurrences. Thus, in using the technique, it may be difficult to detect sudden and previously non-existing changes in data behaviour. This is most obvious in situations of increased randomness and fluctuation in accident occurrences as well as limited availability of historical data. This renders the model incapable of providing satisfactory future transition probabilities in such situations.
This paper develops a grey-fuzzy-Markov industrial accidents forecast model for small or incomplete accident data availability situations using non-differential function grey interval analysis, fuzzy logic, variation conditioning and a state transition approach which aims to capture unique accident occurrence characteristics. The aim of the work is to create a standalone GFM model capable of making accurate industrial accidents forecasts.
The model's development is founded on its ability to recognise accident occurrence variation patterns. These patterns are then decomposed into certain principal pattern components identified in this paper. The results obtained from this knowledge is passed through a fuzzification process and rigorously treated to minimise noise in the fuzzy data. A decomposed state transition approach (DSTA) analogous to the classical Markov state transition approach is subsequently developed and used in detecting future accident vibrations and forecasts are then made in the process.
The validation of the model's existing value and future accident prediction capabilities is done using the in-fitsample and out-of-sample performance evaluation techniques, respectively. It is believed that this novel approach to industrial accident forecasting will aid proper anticipation, planning, control and management of future accident occurrences in industrial organisations on the one hand, and also provide a promising alternative tool to forecasting under uncertain conditions on the other.
The current paper makes a major contribution to the creation of a unique accident occurrence pattern recognition technique based on GFM inferences which acknowledge the significance of uncertainties. As such, the current paper contributes to the discussion on accident uncertainties, which has the interest of accident scholars and also grey-Markov-fuzzy theorists generally.
Industrial accidents forecasting, as argued in this paper, is central to the attainment of industry's stability and a guarantee to survive in the long run since litigation fees resulting from accidents could be reduced to the barest minimum through the adoption of a merit-driven forecasting technique. Nevertheless, the grey-fuzzy-Markov pattern recognition model has rarely been employed to improve forecasting and prediction of industrial accidents in industrial organisations. The authors found a number of papers applying only grey-fuzzy-Markov (Asrari et al. 2012;Geng et al. 2015) in the scientific literature with limited applications to the analysis of electrical and biofuel production, for instance. Industrial accident forecasting has not been tackled in grey-fuzzy-Markov literature. A key issue is that pattern recognition has been under-researched. This shows that the development of grey-fuzzy-Markov pattern recognition framework and the philosophical theory behind it in the context of industrial accidents is a sure gap filled in accident literature.

Methodology
The motivation for the creation of the grey-fuzzy-Markov pattern recognition prediction (GFMAPR) model arose from the observation on preliminary analysis that randomly summative and multiplicative relationships existed between industrial accidents data at different points within an existing data set. The need for the use of fuzzy logic was obvious as it was clearly difficult in employing classical mathematical approaches in understanding such data relationship. 1. Available historical data are randomly occurring and of non-stagnant pattern occurrence feature. 2. Information available for analysis is unique to that system and different in characteristics and behaviour to that of other systems. 3. There is always a summative or multiplicative variation relationship or both existing within any available historical dataset. 4. Second-level variations process has strictly non-static characteristics. 5. A time invariance nature of data exists (Sullivan and Woodward 1994).

GFMAPR: the concept
To be able to develop GFMAPR, two grey-fuzzy-Markov analysis methods, namely summative data relationship (SDR) analysis and multiplicative data relationship (MDR) analysis were carried out on two differently prepared versions of historical data. Grey probable forecasts were subsequently generated from the SDR forecast interval and cross-checked with MDR forecast interval expectations. Based on a set criterion of acceptability, probable forecasts which fell within SDR and MDR interval intersection space were further analysed using a whitenisation procedure to produce the crisp forecast. An outline of the GFMAPR concept is presented in Fig. 1. The procedure for determining the SDR and MDR will be discussed independently in subsequent sections of this paper.

Procedure for the SDR determination
The analysis to determine the SDR is discussed in this section using the outline in Fig. 1.

SDR preparation
Data preparation is the first stage in the SDR process. This stage involves the application of the AGO. At this stage, available historical data x ð1;nÞ fx i : i ¼ 1; 2; 3; . . .; ng were converted into a set of values z i fi ¼ 1; 2; 3; . . .; ng by a cumulative summation of their variations. This is the firstlevel variation analysis.

Creation of summative variation states
Based on results obtained, a set of grey states s k was created to accommodate z i . In creating s k , consideration was given to the dynamic nature of x i evidenced in z i To be able to reflect the current characteristics of data, a position influenced state interval x, assumed uniform for all states was created where x is the parameter which is used to adequately express the relationship between x i and x iþ1 . Thus, the accuracy of GFMAPR is strongly dependent on x.
Obtaining this value requires that two values p and r must be supplied in Eq. (14). p indicates the position characteristic of x i in the data set and r is the set partitioning index. p and r have to be determined. Following preliminary investigation of some industrial accident occurrence data, p was taken in this work as a constant and fixed as r was considered a variable and arbitrarily fixed at an initial value of 4. s k were thus created as follows: The state creation procedure is terminated at k max given that s U k ! maxðz i Þ: To reduce the problem of overestimating the terminal state, s Ã ½sÃ ¼ s k ðk ¼ k max Þ is defined as s L 1 ; s U k ðk ¼ k max Þ Â Ã was adopted as the initial universe of discourse in this work. At this stage r ¼ 4 is not considered as the value that provides the most satisfactory interval width x. The procedure that modifies the universe of discourse by searching for the most satisfactory x was undertaken. This is discussed later.

Fuzzification and reclassification of summative variations
The need to locate z i more precisely within s k necessitated the fuzzification of z i . Due to its simplicity and ease of use, the triangular fuzzy membership function was adopted for the fuzzification procedure in this work. Membership functions for derived s k membership classes are presented in the following equations: where lðz i ; kÞ represents the fuzzy membership values of z i in s k and a k are the midpoints of fuzzy sets s k i ¼ 1; 2; 3; . . .; n; k ¼ 1; 2; 3; . . .; k max ð Þ : The actual state in which z k belongs after fuzzification was subsequently relocated as In-fit-sample forecasts produced by this state classification procedure was then obtained using the following equation: where s ; i is the midpoint value of s ; i . Equation (22) expresses the summative relationship which exists within any given set of historical data related to their variations.
Although the approach establishes that a relationship for a historical data set existing within period n, it cannot be directly employed in forecasting for periods existing outside the historical data window due to the time-invariant data assumption made. The rest of the work is directed at the analysis of the time-invariant model towards obtaining GFMAPR future value forecasts.

Second-level variation analysis and classification of degree of variation
Let variation of z i from its current state s ; i of state number indicates the periodic change in first-level variation within data. Subsequently, the relationship that exists within d i;iþ1 was investigated.

Fuzzy classification of the second-level variations
Preliminary investigation led the classification of the second-level variation (change in the first level) variation d i;iþ1 into four fuzzy classes C b ðb ¼ 1; 2; 3; 4Þ namely: small-level variation (SV), small to medium level variation (SMV), medium to large level variation (MLV) and large level variation (LV). Fuzzy classes based on a trapezoidal membership function d i;iþ1 Z were created for these linguistic classes.
In addition, employing the maximum membership principle, the actual fuzzy class in which d i;iþ1 belongs was obtained as Employing Eq. (29), crisp representatives, Dðd i;iþ1 Þ of C d i;iþ1 À Á were subsequently employed for further analysis.
Expressions (31), (32) and (33) were developed to account for situations of shock occurrences. In such cases, it is possible for intermediate variation classes to be nonexistent, in the presence of higher variation classes. The relations function in smoothening variation levels.
It can be observed that although the SLV was fuzzified, crisp values were employed as the representation of each linguistic class or interval. This was made necessary due to the need to undertake grey analysis on the set of fuzzy classes. which d i;iþ1 ¼ 0 do occur and non-zero D d i;iþ1 À Á equivalents must exist for such situations. To surmount this challenge, smoothing procedures for three different second-level non-variation scenarios were introduced. These are presented in ''Appendix A''.
At the end of the smoothing procedure, let the smoothened values of D j be represented by D These are treated in this section. Note that the term ''pattern swing'' will be used interchangeably with VPCSPs in the course of this discussion.

Identification of data variation principal component pattern swings
Based on the preliminary analysis of several industrial accidents historical data and information, the observation that industrial accidents mostly exhibit randomly trending or fluctuating characteristics or a combination of both were made. Another major feature also observed was that of the presence of various degree randomly occurring shocks within data. Figure 2 shows a cross section of real-time industrial accident occurrences.
Employing these observations, five unique fuzzy VPCPSs which exist in any characteristic variation curve were identified as open (O L ), escalation (E L ), exact-closure (C L ), closure-lag U L ½ and closure-lead V L ð Þ. Letting each D ;; j be the periodic swing magnitude for all periods j the various VPCPSs are defined below: 1. Open pattern swing: This is taken in this work as the next variation swing in the time j given that the previous cumulative variation swing magnitude (CVSM)k jÀ1 is equal to zero. O L can exist at variation curve points occurring immediately after two swings of equal magnitude and opposite poles have offset each other. It can also occur as the sum of swing magnitudes D 00 j and k jÀ1 given that both have opposing poles and the latter is of lesser absolute magnitude.
2. Escalation swing: this is a type of variation swing D 00 j that occurs when k jÀ1 6 ¼ 0 with both having the same polarity.
3. Exact-closure swing: this is expressed as the value of D 00 j with magnitude equal to k jÀ1 but opposite in polarity.
. Closure-lag swing: this is the value of D 00 j with magnitude less than k jÀ1 but opposite in polarity.
. Closure-lead swing: this is the value of D 00 j with magnitude greater than k jÀ1 but opposite in polarity Simultaneously, pattern swing occurrence indicators Variation patterns in this work are generally considered to swing along increasing and decreasing directions. Increasing and decreasing pattern swing directions can be positive or negative at any time depending on current opening pattern swing properties. What is important to note is that while some patterns may swing in a certain direction, others may produce reverse swings. On the basis of this understanding, VPCPSs were then grouped according to the similarity in their ability to swing in the same direction given their presence in data. For example, if the most current swing is in the positive direction and favours an escalating pattern, then if a closing-lag pattern is anticipated as an expected future occurrence, its swing polarity must be in the direction opposite to the escalating pattern and its magnitude must result in a decrease in the current CVSM. Table 1 shows the grouping of pattern swings according to the effect of their swing values on the CVSM. The expected future pattern swing values and directions are dependent on the magnitude and direction of the current CVSM (k cur ), the current cumulative pattern swing(CPS) value P cur and pattern swing impulses r I cur ; r I min and r I max :. k I cur and P cur were, respectively, obtained as k cur values are given primary consideration in the determination of future variations swings. When k cur = 0, then a future open pattern swing is expected. P cur and related impulses are employed for future variation swing determination analysis when k cur 6 ¼ 0.
P cur is used in determining the pre-expected future pattern swing direction for each pattern q i r f À Á (Table 2).
For example, given P cur E L cur À Á , that is, the current CPS being an escalation, if q þ P cur ð Þ is the existing current swing direction, then, the pre-expected future pattern direction for a closure-lag U L f will be the reversed polarity of the former.
P cur can only exist for a single VPCSPs However, a situation can occur in which C L n and V L n may exist within that same period due to their overlapping characteristics. In such situations, a preference to obtain P cur and q P cur ð Þ from V L j is usually made. Pattern swing impulses and related parameters are used in the detecting expected future pattern swings. It is the final stage of the future pattern swing determination. These parameters are obtained from adjusting pattern swing values. The procedures for obtaining them differ from one VPCPS to another. The next section is devoted to discussing this.

Pattern swing magnitude adjustments
A rigorous adjustment process was employed in preparing pattern swing values for future swing estimation. Two adjustment procedures were employed. One procedure was carried out on the basis of CPS impulse and magnitude of occurrence, while the other was undertaken on the basis of the most frequently occurring pattern swings. The two procedures are presented next.
SDR parametric estimation and adjustment based on current and maximum cumulative swings After the split of D 00 J into VPCPSs, parameters related to the duration of swings, recognised as being important for SDR future value analysis were obtained. These parameters are the current CPS impulses for escalation E I cur and closure-lag U I cur ; the maximum and minimum variation pattern swing impulses for escalation E I max ; E I min À Á , exact-closure  In addition, parameters related to actual current and current equivalent cumulative pattern swing magnitudes E ÃL cur ; E hL cur , U ÃL cur and U hL cur were also determined. Before these were achieved E L j and U L j were adjusted to become E ÃL j and U ÃL j . This was based on the understanding that due to the fuzzy nature of the variation properties, there exists the tendency for escalation and closure-lag properties to overlap and potentially occur within respective open and closing swings. These properties were identified and extracted from the respective mother set patterns into their respective subsets. The relations for achieving this are presented in ''Appendix C''. The current and current equivalent cumulative swing magnitudes for the adjusted escalating and closing lag swings were subsequently estimated (''Appendix D'').
Furthermore, there was also the need to update the respective VPCPSs values to reflect their magnitudes in the current period. The update was carried out as follows: With all the estimated parameters and adjustments made, future variation pattern swings were expected to occur based on the set of logical rules presented in Table 3.
Adjustments based on most frequently occurring swing values All pattern swing value adjustments made on the basis of their swing value frequencies ðr / J Þ require a similar procedure. However, determining O / J demands a slightly modified procedure. The steps required for obtaining, and V / f are outlined below followed by the modified form of the procedure developed for obtaining O / j . a. Procedure for obtaining any of Step 2: Step 3: replace all r Lb j values having polarities in reverse of q i r L f and convert data into absolute values.
Step 4: identify the values of r ÃLb j which account for 75% of the set on the basis of the most frequent swing magnitude values. Let the set for which the required data exist be fS75g and b i fi ¼ 1; . . .; i Ã g be the members of fS75g.
Step 5: update r ÃLb j to eliminate swing values that do not belong in fS75g Step 6: determination of future pattern swing. This is the final stage of this procedure. The techniques developed and used for estimating r L j is discussed in the next subsection. b. Procedure for obtaining O / j Notice from Table 3 that when the swing existing in the most current period results in a cumulative swing magnitude of zero, then the expected future swing automatically becomes an open swing. In such a situation, the implication is that there is no certainty on the direction of the next variation in swing. In respect of this, step 1 of the above adjustment procedure loses its relevance. Determination of r L max in step 2 is modified thus, Step 3 to step 6 of the previous procedure is maintained and applied as previously discussed.
Future pattern estimation using decomposed state transition approach (DSTA) The DSTA involves the creation of various transition patterns. If any of the created patterns are dominant within r / j then it is most likely that r / cur will transit from its current state to a future state r L f via the dominant pattern. Four different pattern transition techniques were developed for this purpose. The techniques are (i) same state pattern switch, SSPS a1 ð Þ; (ii) cross state pattern switches CSPS a2 ð Þ; (iii) pattern span measure, PSM a3 ð Þ, and (iv) static dominant patterns, SDP a4 ð Þ. Only one of the first three techniques can be employed during any SDR analysis, while a4 is an inclusive technique adopted for detecting overwhelming static pattern occurrences that cannot be detected by the other three. We proceed to discuss how the methods are developed, the procedure employed in choosing the most appropriate technique for determining r L f and the application of the chosen technique.
Application of the concept of Markov transition in developing the DSTA The four future VPCPS determination techniques were derived by investigating series of swing patterns chains produced by periodic changes in r / j values. Let r / j exist as a Markov chain. Also, let c j;jþ1 be the link chain produced by With respect to the future swing pattern estimation techniques, the following link chains were considered in this study: If investigating for a1 and a2 ; j ¼ 1; 2; 3; . . .; n À 3 f g ð52Þ c j;jþ1 If investigating for a3 and a4 ; j ¼ 1; 2; 3; . . .; n À 2 f g Each periodic link chain investigated revealed a type of swing pattern existing within it subject to certain confirmatory conditions. The frequency of swing pattern occurrence F d within each link was subsequently computed for the four techniques developed. Finally, the probability of VPCPS transition from its most current period to the future period was determined by computing the pattern occurrence strength index I d using the link chain analysis. The developed pattern recognition methods are fuzzified component elements of the Markov transition technique and function by determining the next transition state from the previous (only a single step transition was considered in this work). The methods and procedure for achieving this are outlined next.
Methods and procedure for obtaining the most dominant r / j transition technique Inferences were obtained using r / j , W a and W b 2. Investigate the order of pattern swing occurrence, for a1 and a2 ; j ¼ 1; 2; 3; . . .; n À 3 f g ð55Þ 3. Obtain the frequency of occurrences F d for all techniques d.  (2018) 14:455-489  465 represents the terminal points of the pattern swing chain investigated for which the conditions necessary for s d fhg to exist are fulfilled. These conditions are presented in Table 4. 4. Determine pattern occurrence strength index I d À Á ; a2; and a3 were set in order of preference such that if any of its equivalent I d ! 0:75, then such a technique is considered the most appropriate for determining r L f , while other techniques not yet explored are ignored. Figure 3 is a chart which shows this order of preference Determination of r L f At the end of the future pattern swing value estimation technique determination procedure, given that the desired pattern technique(s) had been detected, the following relations were subsequently employed to obtain r L f .

SSPS
2. CSPS r L f fq;gg¼ 3. SDP q: = a,b; g: = all pattern swing values of different magnitudes belonging to set W q 4. PSM r L f determination by PMS involved some steps. First, the weighted mean ratio of r / j occurrence was obtained.
Employing a3 two possible r L f integer values A and B will exist. That is From these R is the interval bound by min r / Situations in which A; B 2 R were also observed to exist. In such cases, five rules were created and used in obtaining r L f from A and B. The rules and relations for obtaining future pattern values are presented in ''Appendix E''.
1. The analysis just carried out concerns situations of . For the situations involving q À r L f , all r L f ðiÞ obtained was converted to their negative equivalent. 2. Recall that in situations k cur ¼ 0 the predictability of the next direction of swing cannot be out rightly ascertained in this work. To account for this degree of uncertainty, the future pattern swing O f at this stage was analysed in both the forward and backward direction. Thus, the future values for an expected open swing became the union of its future values in the positive and negative directions.
Verification expected future pattern swing values Two activities were adopted in assessing r L f ðiÞ to control and reduce overestimation of SDR forecasts. The activities involve verification to establish maximum pattern swing conformity and conformation to swing time packet limits. This activity involves a comparison of r L f ðiÞ values against r L max . This comparison was necessitated by the need to ensure that r L f ðiÞ does not exceed its maximum allowable swing.
Verification of r^L f ðiÞ for pattern swing time packet limit conformity This verification procedure was undertaken on the outcome of the first stage verification. Here, r^L f ðiÞ values were compared against currently existing variation characteristics to ensure that they did not exceed the existing maximum time packet value limits. To achieve this, there was a need to revert k cur to its time packet values. Let v z be a vector in which the time packet values exist. Then the break-up of k cur into its time packet components v z is expressed as: The time packet value limit verification process is unique to different VPCPS type. ''Appendix F'' presents the relations necessary for verifying r^L f ðiÞ to become r ÃL f ðiÞ. This second stage verification is not applicable to future patterns obtained from expected open pattern swings.

SDR reversal, defuzzification and GFMAPR forecast span generation
At the end of the SDR fuzzification and future pattern determination analysis, r ÃL f ðiÞ were reversed and used in generating GFMAPR forecast intervals. This section undertakes a discussion on the interval generation procedure which was carried out in three phases.
Phase1: SDR reversal The reversal process involved the declassification of the SLV classes into their component fuzzy variation states.
Let Z be a set consisting of the union of r ÃL f ðiÞ r : Phase 2: Determination of fuzzified forecast states and subsequent defuzzification Employing the forecast FLV state numbers k f b were subsequently obtained as: b number of SDR forecast intervals were obtained or generated (as the case may be) by locating the FLV fore- was then obtained from the d f b as: The fuzzy FLV states were then de-fuzzified into actual forecast intervalsx f b using a modified form of Eq. (22), the interval bounds were obtained as: Phase 3: SDR forecast interval generation This phase is an extension of the second phase. It considers two major characteristics of industrial accident occurrences. These are steady trends and fluctuations. GFMAPR was created from two interval variants which give consideration to these characteristics. The two interval variants are the unbound generated interval (UBGInt) and the bound generated interval (BGInt). Equations (77) and (78) are employed for generating the intervals for both variants. Prior to obtaining BGInt some adjustments are made tox f b such that, Notice that BGInt limits forecast strictly to the initial universe of discourse. The use of BGInt for GFMAPR prediction will strongly favour fluctuating situations, while UBGint is created to adapt well to trend occurrences where future values are expected to exist outside a non-preempted universe of discourse. A combination of these two adaptations makes up the GFMAPR.

Procedure for the multiplicative data relationship determination
The SDR future value predictions are intervals values and thus cannot be unitarily employed to obtain point value forecasts. To be able to determine actual accident forecasts, a complementary approach which is the MDR is also developed. The MDR employs a procedure somewhat similar to the SDR but the method differs. The discussion of the procedure is carried out in relation to the outline in Fig. 1.

MDR preparation
The input data for the MDR were obtained from indices resulting from a comparative variation analysis of historical data.
Based on the above, the width of change in variation is expressed as m : U;^; QðUÞ þ Qð^Þ ¼ n À 1 f g s k was then obtained using expression (17) with z i and x replaced by z i and x i ¼ 1; 2; 3; . . .; n À 1 f g , respectively. Note that the boundary adjustment for the SDR initial universe of discuss (expression (18)) is not applicable here. With z i and a k duly replaced in relations (19) and (20), fuzzification and location of z i in s k was also undertaken.
With respect to D j , an exceptional consideration was given to situations where, n ¼ 4 # D j ¼ 1 and # D j ¼ 0 fj 6 ¼ 1g. In such circumstances, the single element in D j with value e j j was adjusted by creating other values between 1 and e j j to exist as elements in D j (expression (86)). Creating this exceptional procedure was necessary because of the availability of a single variation packet having the possibility of a large variation space existing within it. Thus, we endeavoured to glean as much information as possible given such situations by exhaustively investigating the major regions within this variation space.

Determination of MDR second-level variation future values
This is undertaken in two steps, namely (i) investigation of z i to ascertain its most dominant variation characteristics and (ii) determination of future variation values on the basis of z i most dominant trend or fluctuating characteristics. Before proceeding, let b i ¼ d i;iþ1 fi À 1; 2; 3; . . .; n À 2g: b T f is the set of MDR future value forecasts obtained from trend dominant b t . b F f is the set of MDR future value forecasts obtained from fluctuation dominant b t .

Determination of data variation characteristics
Two major characteristics were investigated, namely trend and fluctuations. In estimating the degree of trending g T , and fluctuating g F properties of z i ; the relation below was employed: The quad-point data characteristic trace (Edem et al. 2016) was employed to determine the number of trend chains existing in z i (N T ) and the total number of z i quadpoints (N d ). g T and g F were subsequently employed as b i control parameters for the determination r f fr : Future second-level variation value based on more dominant data variation characteristic The procedure for obtaining b T f ðqÞ differs from that used in obtaining b T f ðq 0 Þ. Both procedures are discussed exclusively.
Procedure for determining b T f ðqÞ If g T ! 0:75, then the following applies: . . .; n À 3g the most current trend chain T w : w ¼ 1; 2; . . .; w Ã f g . This is achieved by obtaining T tðyÞ which are various sets of elements in b i having similar polarities.
The difference between B Ã½LB i and B Ã½UB i is always constant and equal to the variation sway magnitude c ð Þ. v : is the most current pattern swing variation position. Given that g T \0:75, analysis on b i was carried out to obtain B Ã nÀ2 c nÀ2 and v nÀ2 as a prerequisite to determining the future sway direction b Fg f . Equations (100), (101) and   Table 5 show the relations used in obtaining the required requisite parameters (Table 6). jðiÞ Determination of the steadiness in current pattern involved an investigation of b 1 i ¼ n À 2; n À 3; f . . .; n À 6g. If more than eighty percent of b 1 belong to exclusively to either the lower variation class (1, 2) or upper variation class (3, 4) then b 1 was adjudged to be currently steady, else the pattern was considered currently unsteady. If i\5, b 1 is also taken to be currently unsteady. The technique for obtaining B Ã nÀ2 ; c nÀ2 and v nÀ2 involved investigating b 1 from progressively from i ¼ 1 and steadily adjusting and obtaining B Ã i ; c i and m i until the point i ¼ n À 2 was attained. Based on the obtained parameters, b Fg The second method gives stronger consideration to a highly fluctuating pattern in which v nÀ2 does not exhibit huge deviation from c nÀ2 :B Ã nÀ2 was not given consideration in this method. b Fg f 2 values were obtained from the following relations: The desired future forward and backward variation values were subsequently obtained as The intention of carrying out the MDR analysis was to determine a region for which SDR forecast interval could intersect its forecast interval. As a result, F g f point values were considered in addition to all regions before it. Thus, the forward and backward future MDR variation forecasts are expressed as Before proceeding to discuss the defuzzification stage of the SDR, let b f ðg 0 Þ be the set which contains all MDR future variation values b f g 0 ð Þ¼b T f q ð Þ[b f q 0 ð Þ;g 0 ¼1;2;3;...;q 0 Ã n o Defuzzification procedure for the MDR At the end of the MDR future variation determination procedure, a defuzzification procedure involving a twostage process was carried out. The determination of the future SLV values is the first process. This was achieved by identifying various d i;iþ1 values belonging to different D j classes given that D j or À D j exist in b f g 0 ð Þ. Point values of d i;iþ1 for each D j found in b f g ð Þ were then obtained using the weighted average technique.
Recall that an exceptional fuzzification procedure was undertaken for situations involving n ¼ 4. The corresponding first-stage defuzzification relations are The second and final stage of the defuzzification procedure was the determination of the future first-level variation fuzzified state number k f j using d Ã j .
Finally, defuzzified MDR FLV values z f j were obtained from corresponding i f j as: GFMAPR forecasting from whitenisation of SDR and MDR future outcomes GFMAPR forecasting is achieved by a whitenisation process developed in this work. This involves the comparison and adoption of the intersection of grey SDR and MDR future forecast possibilities based on the satisfaction of a fixed intersection criterion. First, SDR future values are reformed to exist in the same orientation as MDR outcomes. This is treated in the first part of this section. The second part discusses the development of the intersection criterion. The final part of the section covers the use of the satisfied and unsatisfied criterion situations for undertaking GFMAPR forecasts.

SDR point future values: creation and reformation
It will be recalled that SDR future possibilities were obtained as a set of grey intervalsx f b b ¼ A comparison of the degree of deviation of M Ãf d from z f j was then investigated.

Development of the forecast acceptability criterion (FAC)
The purpose of determining R jd was to ascertain the closeness of forecasts produced via the SDR and MDR analyses. Thus, M f d values for which corresponding M Ãf d values had very close proximities to those of z f j were adopted into the set of accepted GFMAPR potential forecasts subject to the accuracy of the method used in fixing the acceptable proximity conditions for screening all possible forecasts.
In this work, a proximity score method was employed in establishing the FAC. Equation (120) infers that a perfect fit between SDR and MDR forecasts will produce zero width of deviation. An initial score was chosen for this forecast scenario. The scoring method for other scenarios was developed from the understanding that z f j is also a fuzzy value even though its parent state has not been established. z f j was assumed to exist as the midpoint value of a pseudo-state. The regions bounding the pseudo-state were taken as being equal to the state width value x . Any R jd value within À2x R jd 2x was awarded a score e jd within the range 0 e jd 10. Subsequently, it was expected that any R jd existing within À2gx ; 2gx ð Þg [ 1 f g will possess score values in the range 10\e jd 1: Employing these assumptions and adoptions, a general equation obtained from the interpolation of the identified points and equivalent fixed scores was developed for computing e jd for all R jd : Notice that although x was our assumed state width from z f j 2x was employed instead. This was necessary to account for uncertainties and model inadequacies.
For the least deviating forecast to accrue the highest score and the most deviating forecast the least, the outcome from Eq. (121) was employed such that minðR jd Þ was given the score of maxðe jd Þ while minðe jd Þ was awarded maxðR jd Þ. The entire score for R jd was then recomputed.

GFMAPR single value forecast analysis
The analysis to obtain x f nþ1 ðx f nþ1 ¼ x f ðx 1;n ÞÞ is the final stage of the GFMAPR forecast analysis. Two methods used exclusive of each other for obtaining x f nþ1 were employed. The first and more prioritised method considered a situation in which the FAC analysis revealed that 9R jd : À2x R jd 2x . The second method applied when the first condition was violated, that is, 8R jd : ðR jd \ À 2x Þ È ðR jd [ 2x Þ. A procedure which shows how x f nþ1 was determined using the methods is outlined below.

Procedure for GFMAPR forecast determination
i. Set initial value of jðj ¼ 1Þ ii. Obtain R jd fd ¼ 1; 2; 3; . . .; d Ã g. iii. Apply the FAC criterion. Determine max e jd À Á . Also, obtain corresponding M f d max e jd À Á À Á . Test to see if 9R jd : À2x R jd 2x . If this condition does not exist go to vi. iv. Compute the proximity score index C jd for all e jd C jd ¼ 100 Â max e jd À Á À e jd À Á min max e jd À Á ; e jd À Á h i v. If C jd 10, accept M f d into the set of GFMAPR potential forecast X.
Let corresponding e jd also be an element in Y vi. If j\j Ã increase j by a unit value and return to ii. vii. If j ¼ j Ã then, a. If b. 9R jd : À2x R jd 2x , obtain x f as, viii. End procedure It is worth noting that: Search procedure for improving GFMAPR forecast A preliminary investigation was undertaken to observe the accuracies of GFMAPR forecasts with respect to varying SDR sets. To this end, the set partitioning index r in Eq. (14) was replaced with r i ¼ fi ¼ 1; 2; 3; . . .; 1g; where The investigation led to the observation that GFMAPR forecast accuracy had the potential to improve at different r values other than r ¼ 4. Further experimental investigation of model performances on application to different industrial fire accident data evaluated on the basis of the MAPE showed that better GFMAPR performances could be obtained from the number of created SDR sets related to a r ¼ r ð128Þ x Ã is the value of x for which one of the set partition indices r Ã , obtained from an initial search using r i ¼ fi ¼ 1; 2; 3; . . .; qg produces the best GFMAPR performance, measured using some performance evaluation approach. a r is a multiplier for r Ã employed for an extended search for better GFMAPR performances within regions covered by multiples of r Ã .
As an example of how r Ã was determined in this work, GFMAPR out-of-sample forecasts were carried out at various set partition index values using three industrial accidents occurrences from three literature sources. This was achieved by obtaining the best model MAPE performance at r i ¼ fi ¼ 1; 2; 3; . . .; 5g. r i for which model best forecast was obtained was recorded asr. Subsequently, the search for improved model performances was widened by further applying the model to obtain forecasts for set partition index values of vr AE i fi ¼ 1; 2; . . .; b r À 1g. v is an integer which represents multiples ofr v ! 1 f g. r Ã ¼ vr. Figure 4 shows the points of troughs p which indicate model best forecast performances at set partition index regions r Ã þ i ¼ f1 i\r i g (Table 7).
Two observations were made from Fig. 4. The first is that although the regions where GFMAPR produced relatively high-accurate forecast did not exist as r Ã multiples, nonetheless a search around the region of r Ã a r showed potential to improve GFMAPR performances. Second, it was also observed from Type 2 and Type 3 data analysis that the strict use of r ¼ 4 did not always guarantee the best forecast performance of GFMAPR. The need to improve model performances for all types of industrial accidents data based on these observations justified our interest in developing the GFMAPR search (S-GFMAPR) model. Subsequently, the GFMAPR variant in which r ¼ 4 is referred to as the GFMAPR non-search (NS-GFMAPR) model.

Development of the comparative performance index for detecting best S-GFMAPR forecasts
To obtain S-GFMAPR forecasts that are not strictly hinged on a single performance evaluation measure, a comparative performance index (CPI) was developed. The CPI is developed by combining GFMAPR out-of-sample results. This result was determined by the use of a single horizon Fig. 4 A plot of GFMAPR (MAPE) performances for a set partition index range of 1-25 using data presented in Table 6 Table 7 Description of points of high GFMAPR performances of Fig. 4 in terms of the first trough point r Ã Data Type 1 P1: r Ã Ind Eng Int (2018) 14:455-489 473 recursive rolling forecast approach (see Sect. 8.2) which was used in computing model MAPE, MAE and MSE result determined during the search process. It was used to distinguish between S-GFMAPR performances computed using two most closely following set partition indices. The superior S-GFMAPR performance was taken as that which produced a lower CPI when any two sets of performances determined within a search region (r) were compared with each other. For each round of comparison t r ft r ¼ 1; 2; 3; 4; . . .; t Ã r g the superior performance evaluation index CPI Ã t r ð Þ was obtained as MAPEðbÞ MAEðbÞ MSEðbÞ ðb :¼ t r À 1; t r Þ t r : is the current period within r region in which GFMAPR results have been evaluated. t r À 1 : is the previous period within r in which GFMAPR results have been evaluated.
N t r ð Þ : is the number of periods/rounds available for search investigation.
CPI Ã t r ð Þ : is the lower CPI value obtained after comparison in the period t r .
As a result of the limitation placed on the minimum data size requirement for GFMAPR, the out-of-sample performance index which can be obtained from Eq. (131) can only be computed if an available number of historical data is greater than four. In situations where the historical data number available is equal to the minimum, the S-GFMAPR approach cannot be employed. The NS-GFMAPR becomes useful for forecasting in such situations.

Algorithm for undertaking S-GFMAPR forecasting
Based on conclusions drawn from our preliminary investigation, a procedure for carrying out the search GFMAPR forecasting was developed. The procedure is outlined below using the following steps.
Also, obtain corresponding x f nþ1 CPI Ã r ð Þ Step 6: Update CPI ÃÃ ðrÞ and obtain corresponding Step 7: Step 8: If PTS ¼ 4 then GFMAPR forecast x f nþ1 at the end of the search procedure is r Ã a r ðt r Þ r Ã a r r Ã a r þ 1 r Ã a r þ 2 r Ã a r À 1 r Ã a r À 2 Otherwise, if PTS\4 return to step 4 and continue the procedure.

Results and discussion
The GFMAPR forecasting procedure could be cumbersome and time consuming to undertake manually. In this regard, a Visual Basic.Net program was created to execute the multiple procedures contained in the model.
Model performances were investigated on two fronts, namely accuracy of the developed set creation and partitioning method which was measured using the insample (trained model) prediction evaluation, and the forecast capability of the model which was investigated using the out-of-sample forecast method.

GFMAPR data set partitioning accuracy in model training
In-sample fitted results obtained using the GFMAPR (SDR) set partitioning technique was compared using two industrial accident data (Zheng and Liu 2009;Kher and Yerpude 2015) and one traffic accident data (Arutchelvan et al. 2010) The traffic accident data were adopted on the  basis that traffic accident occurrences also exhibit random and uncertain patterns. In addition, the traffic accident data have been employed in the analysis of several set partitioning models (Lee et al. 2007;Jilani and Burney 2008;Egrioglu 2012;Kamal and Gihan 2013). Table 9 shows actual in-sample (trained model) predictions for one of the accident occurrence data that were analysed, while Tables 10, 11 and 12 show the NS-GFMAPR and S-GFMAPR prediction fitness to data used in building the model in comparison with results obtained by the established models previously mentioned. It can be observed from the presented results that the set partitioning technique developed in the GFMAPR model produces the best fit to data when compared to various results obtained from commonly employed models employed for accidents prediction.

Forecast capability of GFMAPR
It can be observed from Table 9 that the apart from predicting trained data outputs, the developed model can also be employed for making future predictions (referred to in this paper as forecasts). However, the model is limited to making forecasts for single horizon (one step ahead) only. That is, given an available dataset of n fn ! 4g occurrences, GFMAPR can only forecast accident occurrences for the period n þ 1 and not beyond that. Thus, the model may not be suitable for making multiple horizons forecasting or model out-of sample tests except in situations where bootstrapping methods are utilised. As a fallout of this limitation, splitting a dataset into model training and testing portions become unnecessary. Thus, for a given dataset of n occurrences, n number of data points are also required for training GFMAPR. This property of the model imposes a restriction on proper model validation with respect to the out-of-sample prediction or forecast capability of GFMAPR. In overcoming these shortcomings, the out-ofsample forecast evaluation method was adopted for model evaluation and validation.
The out-of-sample forecast evaluation exploits the single forecast horizon capability of GFMAPR using a recursive rolling mechanism. This means that for any data sample of size n with occurrences in horizon k xðkÞ k ¼ 1; 2; 3; . . .; n f g ½ , the forecast x f ðkÞ was obtained for a trained split data set x ð1;kÞ ðk ¼ 4Þ and tested against x ð1;kÞ . The split data set was subsequently increased by a single data point and the training and single horizon forecast (testing) process repeated for k ¼ 5; 6; 7; . . .; k ¼ n À 1.
The PE indices ðPE mij Þ for GFMAPR m ð Þ in terms of the absolute error (AE), absolute percentage error (APE) and square error (SE) were determined using x f ðiÞ and x iþ1 fi ¼ k; k þ 1; k þ 2; . . .; n À 1 : k ! 4g. PE oos mij for each model variant was then obtained as, The capability of the GFMAPR variants for single horizon forecasting was subsequently investigated by comparing their forecast performances with results produced by six forecasting models with forecasting capability commonly used or designed for use in industrial accidents forecasting. The established models employed for the analysis were namely: ARIMA, ESM, a three-point M.A model, GM (1,1), Grey-Markov model (Mao and Sun 2011) and DPEWTA (Edem et al. 2016).
The model validation and evaluation involve a comparison of the results of the GFMAPR model variants with those of the mentioned models on their application to some industrial accidents data using the out-of-sample forecasting approach previously discussed.
In carrying the evaluation, a PE weight indexm was created by converting the MAPE, MSE and MAE values obtained for all the models into relative weight values C mj (Eq. 141). The sum of these relative weights was then deployed as the PE index for the model variants as well as the compared models.  NðMÞ is the number of models employed for comparative analysis} NðJÞ is the number of performance evaluation methods employed for comparative analysis} This approach was necessary as a result of the need to utilise a single PE index which exhibited the combined characteristics of the MAPE, MSE and MAE in the evaluation of the models.
Using the out-of-sample forecast evaluation approach, the forecasts made by the developed and compared models when applied on the three accidents historical data sets previously employed in Sect. 8.1 were obtained. The out-of sample forecasts for the three industrial accidents data employed in this work are as shown in Tables 13, 15 and 17. The forecast performance evaluation results obtained from the three accidents historical data sets using this technique are shown in Tables 14, 16 and 18.  -GFMAPR  S-GFMAPR   1990  91  --------1991  80  --------1992  107  --------1993  101  --------1994  93  94  97  96  100  118  124  94  94   1995  91  95  93  95  97  103  106  101   It was observed from the results presented in Tables 14, 16 and 18 that in terms of them which combines the accuracies of the models MAE, MSE and MAPE into a single evaluation score (Eq. 135); S-GFMAPR produced the best out-of-sample forecast results in all but one of the data set analysed. In terms of the performance score,   -GFMAPR  S-GFMAPR   1974  1574  --------1975  1460  --------1976  1536  --------1977  1597  --------S-GFMAPR produced the highest score in two of the three analysed data with a minimum and maximum positive relative score difference of about ten percent and seventy percent, respectively. The model, however, produced inferior results when compared with grey and Grey-Markov models on application to the chemical plant accidents data.
The general results obtained indicate that although S-GFMAPR may not be suitable for forecasting all forms of industrial accident occurrences, it does exhibit superiority over other compared models when applied on data with characteristics suitable for S-GFMAPR forecasting.
However, the computation time (CT) for GFMAPR models using a 2.13 GHz Intel Pentium P6200 CPU with 4 GB of RAM is generally higher compared to all other models. Nonetheless, the CT range of 1 and 2.5 s can be considered as acceptable tradeoffs for the models excellent forecasting capabilities. As expected, NS-GFMAPR produced results less superior to those of S-GFMAPR. However, when compared in terms of all eight analysed models, its general performance score point ranking fell between 3 and 4. This implies that although it is less effective than S-GFMAPR, it can still be employed for forecasting. NS-GFMAPR will be more suitable in situations where time efficient forecasting is a necessity, and when available historical data points are limited to four. To further ensure the validity of the GFMAPR variants, the models were applied on five more real life data related to industrial and traffic accident obtained from literature sources (''Appendix G''). The results obtained further strengthen the validity of NS-GFMAPR and S-GFMAPR as   S-GFMAPR   2000  350  --------2001  347  --------2002  437  --------2003  260  --------suitable models for forecasting fluctuating, random and uncertain occurrences.
Investigation to ascertain GFMAPR forecast capability given varying historical data sizes and shocks Although preliminary investigation showed that the GFMAPR model possessed huge potential for more accurate small data sample size forecasting, a more exhaustive investigation was undertaken to ascertain this observation. Consideration was given to small (4-15), medium (16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35) and large ([ 35) sample size of data. Due to the unavailability of industrial accident occurrence records of very large sizes, simulated data were employed. Thirty data sets, each made up of seventy periods was simulated. To ensure that the simulated data possessed real industrial accident occurrence characteristics, properties of historical industrial accidents data obtained from the literature were investigated and employed for the simulation. The properties investigated for include the degree of shocks and the frequency of fluctuation in occurrences. Each data set characteristic was simulated using the information turbulence index ðhÞ and the Quad point characteristic trace QPCT approach (Edem et al. 2016) such that each possessed some degree of shocks within the range of 0\h 0:4 and fluctuation index ðg F Þ ranging between 0\g F 1 (Eq. 89). The recursive rolling single horizon out-of-sample forecasts were undertaken for S-GFMAPR, ARIMA, ESM, MA, Grey-Markov and DPEWTA models.
GFMAPR forecast MAPE performance relative to each of the models employed for comparative analysis R mdk ð Þ was then computed as A graphical plot of R mdk against k was then undertaken for each m.
At the end of the analysis, it was observed that GFMAPR exhibited forecast superiority over all compared models with respect to small data size forecasting. It also showed superior performance over grey, grey-Markov and DPEWTA models with respect to medium and large data sizes. However, ARIMA and ESM exhibited superior forecast performance over GFMAPR within data size range between 28 and 70 (Fig. 5). Note that negative values of R mdk imply superior GFMAPR performance.
Furthermore, the simulated data sets were separated into five equal groups based on their degree of variation characteristics. Their MAPE performances were also observed as they changed with data size. The MAPE for GFMAPR relative to each compared model showed that GFMAPR forecasts increasingly improved with increasing degree of variation in data. For example, Fig. 6 shows the relative GFMAPR-ARIMA forecasts for increasing h values. It can be easily seen from the figures that at values of h [ 0:3, GFMAPR models become more reliable for undertaking forecasts. A summary of the results obtained from this and analysis is presented in Table 21.
The mean ð FÞ and standard deviations ðr F Þ of the GFMAPR forecast accuracies on its application to the simulated data samples given different variations and size of data samples are presented in Table 20. It is believed that these results are statistically valid since thirty datasets were simulated for each data point that makes up a sample size class. The obtained results further strengthen the capability of the model for accident forecasting.
It should be noted here that the all forecast results obtained are acceptable within a single forecast horizon which is the focus of this work. The validation of the models with respect to their capability to accurately forecast industrial accidents was also investigated. This was achieved by the application of the model to a case of real-time industrial accidents occurrence data.
The historical information for the case study represents the accidents data of a cocoa processing firm operating in the south-west region of Nigeria. The data obtained covered operating periods within 2003-2015. The required data for the model were collected from the vetting of accident records of various departments operating in the organisation. During the process of data collection, the accidents were observed to occur from sources such as actions and inactions related to poor housekeeping, lack of proper communication during equipment repair and maintenance and delay in the elimination of identified hazards.
GFMAPR forecasts results obtained from the analysis of the case study are presented in Table 19. Based on Lewis's subjective MAPE interpretations (Ofori et al. 2012), the model variants produce MAPE results that indicate reasonable forecasting. In addition, considering the data size (14 periods) and the degree variation in occurrences (h ¼ 0:2592) the GFMAPR forecasts evaluated in terms of the MAPE can be considered satisfactory since it was found to lay within an expected range of F AE 1r F ð4 k 15; 0:2 h 0:3Þ (See Table 20).
The MAPE result may be deceptive due to its scale sensitivity in the presence of a low volume of accident occurrences which is the case of the organisation under study. In line with Stellwagen (2011) recommendation on the tests of model validity in such situations, the model forecast evaluated using the MAE shows forecast accuracy that lies within the accident occurrence range of AE3 and AE2 for NS-GFMAPR and S-GFMAPR, respectively. This  result further strengthens the accuracy of the model (Table 21).

Advantages and disadvantages of GFMAPR
The GFMAPR model proposed in the current paper has a number of advantages put side-by-side of earlier developed models. Here, the GFMAPR is noted as an excellent model in the absence of exhaustive data. In fact, it generally exhibits superior short and medium data length forecasting quality over all models employed for comparative analysis. This is substantial as a merit, and one of the strongest attributes of the model. The GFMAPR does not require many technicalities in execution. Interestingly, GFMAPR requires only the input information. Data pattern analysis and detection are executed by the model. Another merit of the model is that it does not require parametric input to describe data characteristics. It is also noteworthy that GFMAPR model exhibits improved forecasting ability overall compared models when applied on data characterised by strong variations. Regrettably, even with the enormous advantages of using the GFMAPR model over other competing ones, it has a number of weaknesses. First, the model is incapable of undertaking out-of-sample forecasting for periods longer than one forecast horizon. Second, computational experience during the analysis of the model results showed that it is difficult to execute manually without the aid of a computer device. A third drawback of the model is that it is less efficient in terms of computation time when compared with other established models analysed in this work. It was also observed during the model analysis that GFMAPR model was less effective when employed for large size data forecasting. This is a disadvantage. Nevertheless, this disadvantage may be ignored since the primary aim of developing the model was hinged on the need to make accurate forecasts in the absence of extensive data sizes.

Conclusions and future works
A fuzzy set classification based method which employs grey, Markov and pattern recognition concepts for forecasting industrial accident occurrences has been developed. The in-fit-sample and out-of-sample forecast performance of the model was investigated in comparison with some forecasting models frequently employed in the industry by their application on three real-time industrial accidents related data. In addition, the forecast capability of the model with regards to data sizes, variation and fluctuating characteristics were also investigated with the use of simulated industrial accident occurrence scenarios.
In this communication, our purpose is to present an understanding of the first analysis from both the theoretical perspective and a practical application viewpoint of the conceived theory of grey-fuzzy-Markov pattern recognition model for accident forecasting under reasonable assumptions. This investigation broadens the horizon on accident forecasting in industrial settings by showing that accident data could be read into patterns and successfully decomposed and then using the transitional Markov property to detect the possible vibrational direction of future accidents. The GFMAPR model was developed to enquire about the linked characteristics of grey, fuzzy, Markov and pattern recognition.
Through several case examinations against published literature a real-time data and, the main conclusions are as follows: • GFMAPR shows excellent forecast capability for data characterised by fluctuations, variations and randomness; • GFMAPR is quite effective for forecasting in the presence of small and medium sized data availability; • GFMAPR is able to make relatively high accurate forecasts when compared against models frequently employment in the industry for forecasting; and • GFMAPR is thus suitable and reliable for industrial accident prediction.
This paper has made important contributions to research as well as practice. On the first note, an original method based on grey-fuzzy-Markov pattern recognition model was developed to track uncertainties, imprecision and randomness in historical data. Secondly, a new summative data analysis, which involves data preparation, the creation of summative variation data relationship (SDR), fuzzification and reclassifications of data based on the degree of variation was developed. A practical step in evaluating the association among data elements was established for the first time. Similarly, the multiplicative data relationship (MDR) analysis was established as a first-time contribution in accident forecasting literature. This investigation, in search for a procedure in enhancing GFMAPR forecast, developed a comparative performance index for detecting the best S-GFMAPR forecasts for the first time in literature. Furthermore, it is the first report in literature, on accident forecasting, to have an algorithm for understanding the novel S-GFMAPR forecasting, and this has been clearly detailed out in the current paper. Now, considering the contributions of this paper, it hoped that this report will serve as a reference document, as procedural paper that guides further investigations in further probing the details of the foundations and in so doing; we will have a rich number of papers that may develop into an area in industrial accident forecast. Consequently, this study has a practical importance for industrial systems to work out and execute an accident forecasting model that accounts for uncertainties, imprecision and randomness in their systems. Based on this, industrial organisations would acquire the will of monitoring and planning fully for any accidents that would occur in their systems.
Researching to unveil the potential of grey-fuzzy-Markov models with respect to accident forecasting is definitely a new topic in the safety arena. There are lots of opportunities for future research on the topic. Future investigations could be directed at improving the industrial accidents forecast by studying how combined or hybrid models can be developed using GFMAPR for the purpose of producing more improved forecast performance. In addition, further research could be made towards the use of GFMAPR for multiple horizon forecasting. It may also be worthwhile in finding out the performance attribute of realtime data collected from rarely studied systems on accidents such as maintenance functions as much as accidents appear to be due to poor maintenance actions in organisations probably due to inefficiency, poorly trained skilled workers some other causes.
In this communication, theoretical results have been provided based on the attributes of the GMFAPR model using the combined characteristics of grey, fuzzy, Markov and pattern recognition. Forecast performance of the GFMAPR model was compared model for the varying degrees of data variation range as 0 h 0:5, 0:05 h 0:1,0:1 h 0:2, 0:2 h 0:3 and 0:3 h 0:4 while the data size range according to small 4 k 15 medium 16 k 35 and large 36 k 70. It was shown that forecast results using GFMAPR were progressively enhanced with small and medium size data and with elevations in the degree of fluctuation and variation attributes for the explored data. The theoretical model outcomes were demonstrated through scrutinising the performance of models as they were related in comparison with GFMAPR using models as ARIMA, ESM, MA, MSGM and DPEWTA. Thus, from the quantitative analysis using data sets, the forecast with support from GFMAPR are reasonable with the aforementioned alternatives presently in vogue.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Let
D d i;iþ1 À Á exist as row vector D j fj ¼ i : i ¼ 1; 2; 3; . . .; n À 1g. Also, let the point for which D j ¼ 0 be m. In addition, set the initial smoothing cycle counter t ¼ 0 such that the vector within any procedure cycle is expressed as D j ðtÞ. The smoothing procedure can be carried out as follows: Return to 1 ii. If m 6 ¼ 1 and m 6 ¼ n À 1 then If 9 D j ðt À 1Þ : D j ðt À 1Þ [ 1 Obtain D q ðt À 1Þ and D w ðt À 1Þ Then, Return to 1 iii. If m 6 ¼ 1 and m 6 ¼ n À 1, then if max . . .; g Ã ; D g Ã 6 ¼ 0 D m t ð Þ. Return to 1.
b. If m does not exist, end procedure.
It is worth noting from the smoothing procedure that one or more of these static variation scenarios can exist in D j given single or multiple non-variation points. The approach here is to employ the procedure in eliminating each nonvarying point. Thus, for T non-varying points to be eliminated, then T cycles of the smoothing procedure will be applied. g Ã is the termination point of the procedure for one cycle.

Appendix B
Relations for obtaining current, maximum, and minimum variation pattern swing impulses i. Current pattern swing impulse: ii. Maximum, and minimum pattern swing impulses for escalation and closure-lag: r I min ¼ minðr i Þ:r I max ¼ maxðr i Þ; ðB4Þ iii. Maximum, and minimum pattern swing impulses for exact-closure: :j ¼ 1; 2; 3; . . .; n À 1 Appendix C Adjustment relations for escalation and closure-lag patterns to account for overlap properties where, Table 22 Values of r hL cur adopted for various r hLUB T q þ C L j ; F C L j ¼ 1; C L j C L : > : : Appendix D Procedure for the estimation of current and equivalent cumulative swing magnitudes for escalating and closinglag patterns (see Table 22 ii. Determination of E hL cur andU hL cur : The following steps are employed in estimating E hL cur andU hL cur is outlined below using the following steps, Step 1: Step 2: Using time weights, adjust absolute r 0 b to become Step 3: Obtain H and h, respectively, Step 4: Obtain r 0w b variation index: r hL cur will be a value existing between (r hLUB cur and r hLLB cur ): Note that the elimination approaches are arranged in their order of preference of use. If a more preferred method is explored and found to be suitable for use, then the other methods are subsequently ignored. Otherwise the next preferred method is explored.