Wet precipitation (WP) is a diffusion-controlled synthesis method, which is often used for synthesizing such compounds as hydroxyapatite (HAp). Since the process is limited by diffusion, the choice of a diffusion model becomes a critical aspect. In this simplistic assessment for a preliminary evaluation of the diffusion model applicability, the Ginstling-Brounstein (GB) equation is chosen and analyzed for the case of spherical particles. The nominal kinetic constant K is a parameter in GB model which describes diffusion and is related to the effective molecular diffusivity. When the value of K is known, it becomes possible to predict the required time to achieve desired conversion and design the synthesis accordingly. The GB model is assessed mathematically using simulations, a parametric study and Yates analysis (2n factorial design). Parameters chosen for a preliminary study are in the range of characteristic values for a laboratory-scale WP synthesis of HAp and are thus representative for the application of the model to practice. It should be noted that the analysis is simplistic and is meant to provide only preliminary information for future research, requiring experimental validation.
The Ginstling-Brounstein (GB) model is applicable for describing diffusion-controlled processes, and its practical use has been investigated and reported for a wide range of applications [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18]. Known applications that mention the GB model for describing the kinetics include the hydration of the portland cement particles , solid-state synthesis of lanthanum manganite controlled by the three-dimensional solid-ionic diffusion , solid-state synthesis of the compound Zn2.5VMoO8 , formation of magnesium ferrites , dehydration of the iron(III) phosphate dihydrate [5, 6], thermal degradation of carbohydrate polymers , oxidation and decomposition of three-dimensional braided carbon fibers , alkaline hydrolytic decomposition of the uncolored chromium leather wastes , thermal stability evaluation of silver sulfathiazole-epoxy resin networks , oxidation of boron powder , wollastonite fibre dissolution in acetic acid aqueous solution , xylan pyrolysis , thermal degradation of DGEBA epoxy crosslinked with natural hydroxy acids , and pyrolysis mechanisms of lignin-PVA blends . The model has also been applied to some extent in the study of thermal modes of heterogeneous exothermic reactions  and the study of inclusion complexes in supramolecular host-guest architectures , as well as a methodological study on determination of kinetics from DTA and TG curves .
Wet precipitation (WP) is one of the most widely used hydroxyapatite (HAp) synthesis methods and is based on the reaction of calcium hydroxide with the orthophosphoric acid diffusing into the particle. Hydroxyapatite is the main mineral constituent of bones and teeth, and has found its application in various fields such as tissue engineering, orthopedics, prosthetics, drug transport and environmental remediation [19, 20, 21, 22, 23]. The orthophosphoric acid is added dropwise into the suspension. The reaction can be written as in the reaction Equation 1.
In this work, the diffusion model is analyzed in the context of diffusion of orthophosphoric acid into the calcium hydroxide particle in a water suspension. However, it has to be noted, that in the studied case, the chemical reaction is much faster than the diffusion, thus, the orthophosphoric acid is diffusing into the reaction product (hydroxyapatite)  in reality. The model is somewhat simplified, since it does not account for the presence of byproducts which are normally created in the case of hydroxyapatite synthesis .
Another important assumption is that particles can be considered as spherical, which is often only an approximation. The coefficient D is the effective diffusivity. It describes the diffusion in the porous material, and is macroscopic in its nature, e.g., describes not the individual pores, but all the available pores in total .
where τ is tortuosity; ϕ is available porosity; σ is a contraction factor.
The available porosity is equal to the total porosity excluding pores, which are not accessible to the diffusant due to geometrical dimension restrictions, and excluding pores, which are not connected with the rest of the pore system and thus are closed. The contraction factor describes the decrease in the diffusion rate due to the increase of viscosity in the close proximity of pore walls (Renkin’s effect) [28, 29].
Since the diffusant reacts and diffuses into the particle, the resulting differential equation in the spherical coordinates in the steady state conditions can be written as in Equation 3 . The diffusant’s concentration at the particle’s surface (r = 0) is CAS.
where n is the reaction order; kn is the reaction rate constant.
Based on the Equation 4, it can be deduced that the solution structure does not allow to separately obtain kn and D, simultaneously, but rather only the ratio. Taking this into account the simplest model is chosen and is further described. In this work, GB diffusion model is chosen, which is based on earlier Jander’s work [30, 31].
Diffusion in the spherical particles is a limiting stage in various chemical processes, including the heterogeneous reaction of orthophosphoric acid with the calcium hydroxide. Describing such processes in a general sense, a compound A reacts with a compound B, resulting in a compound AB, whose layer thickness continuously grows during the process. Furthermore, the compound A diffuses into the compound’s B surface through the layer of compound AB with such rate that is much lesser than the chemical reaction rate of the A with the B .
Since the outer resistance to diffusion is significantly less than the inner resistance of AB, compound’s A concentration in the plane dividing A and AB can be considered constant. In the plane dividing AB and B, the concentration of compound A is always equal to zero due to much higher chemical reaction rate between A and B than that of diffusion .
As the starting point in GB model a classical diffusion equation for spherical particles is used, represented by Equation 5.
where the is the proportionality coefficient when concentration is in the molar units; n is the stoichiometric reaction coefficient which represents compound’s A quantity in moles reacting with a single mole of compound’s B; μ is compound’s A molecular mass.
As one of the solutions intermediate forms to this model the Equation 10 is obtained.
Equation 10 describes the relationship between the thickness of the reacted layer x and the elapsed time τ. It can be deduced that the rate of the change in layer thickness of compound’s AB continuously decreases in time at values from 0 up to 0.5, but then symmetrically increases at values from 0.5 up to 1.
Since it is extremely challenging or even close to impossible to experimentally determine the thickness of the reacted layer in time, Brounstein and Ginstling proposed to use the conversion G in place of x, where G varies from 0 to 1.
In the G – τ coordinates, model can be represented by the Equation 11.
where K is the nominal kinetic constant, and is related to the effective diffusivity as represented by Equation 12.
where Φ(y) is the error integral defined as in Equation 13.
Equation 12 can be solved numerically and the value of D can be obtained. However, in this work, from the practical standpoint, it is more efficient to look into the ways of experimentally evaluating the nominal kinetic constant, since effective diffusivity includes both the diffusivity and reaction kinetic terms, separate evaluation of which is an unnecessary complication in this case and thus to be avoided.
The final GB model is given by Equation 14.
where G is conversion (-); K is nominal kinetic constant (m2/s); τ is time (s); R is particle radius (m).
The nominal kinetic constant K is a parameter in GB model which describes diffusion and is related to the effective molecular diffusivity. When the value of K is known, it becomes possible to predict the required time to achieve the desired conversion and design the synthesis accordingly.
Experiment simulations and analysis methods. The heterogenous reaction occurs at the Ca(OH)2 particle surface. The diffusion model is assessed using the aforementioned GB model. The mathematical model is assessed using experiment simulations. The full output generated during the simulations is available in the Appendix of this work. The particle radius is assumed to be in the range of 0.5 – 1.5 μm. It is known, that in practice it is common for the wet precipitation (WP) synthesis of hydroxyapatite to long for about 5 to 7 hours. It is then reasonable to assume that at this elapsed time the reaction is close to a completion and the conversion can be assumed to be in the range of 0.90 – 0.95. Using the GB model, it is then possible to estimate the value of K.
The experiments are simulated, and the effect of K and R on conversion G is analyzed, as well as the effect of parameters G, R, τ on the value of K in calculations, with the help of statistical methods. 45 simulations were performed in total and are reported in the Appendix of this work. The data is analyzed using 2n factorial design using Yates Algorithm, also known as Yates Analysis, in order to identify the percentile effect of parameters on the value of K in the model. The model in the form of Equation 12 is used.
Ethical approval: The conducted research is not related to either human or animals use.
3 Results and discussion
3.1 Experimental simulations
Numerical simulations of 18 experiments were performed, that is for the all possible combinations of extreme cases and also a few medium cases for particle radius and elapsed time within the boundaries of the first assumption as represented by the Equation 15. Values chosen for a preliminary study are in characteristic range for a laboratory-scale WP synthesis of HAp.
The algorithm of the simulations is as follows: the input parameters (particle radius, conversion and elapsed time) are inserted in the model. The value of K is obtained. Knowing the value of K, and setting G to a value from 0 up to 1 with a step of 0.05, it is then possible to calculate the time required for obtaining the respective conversion.
The results of these 18 simulations (S1 – S18) can be found in the Appendix of this work. The example of the simulation output with the input parameters being R = 1.0 μm, G = 0.9, τ = 7 h is shown in Table 1.
|1 - (2/3)G - (1-G)^(2/3) = K*τ/R^2|
|G||τ, S||τ, h|
Obtained values of the nominal kinetic constant K are summarized in Figure 1.
As shown in the Figure 1, it can be deduced that the expected value of K is in the range of about 1∙10-19 to 3∙10-17 m2/s. It should also be noted, that the larger the R, the larger the scatter of possible K values. Based on the 18 simulations, the results are analyzed via descriptive statistics and the average value of K in the expected range is reported in Table 2. Based on the aforementioned assumption (Equation 15), the expected value of K in the practical conditions of wet precipitation synthesis of hydroxyapatite is 1.14∙10-17 ± 4.33∙10-18 m2/s (95% probability), which results in 37.87% uncertainty.
3.1.1 Influence of the particle size and nominal kinetic constant on the conversion
Analyzed particle radius influence on the conversion change in time, using the obtained average value of K = 1.14∙10-17 m2/s from the aforementioned 18 simulations (S1 – S18). Additional 6 simulations (R1 – R6) are performed, where the value of K is fixed, and 6 different R values are chosen: 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 μm. Note, that radii higher than 1.5 μm are also used here, in order to obtain more data points and thus a fuller picture of the particle radius influence. Elapsed time in order to reach a specific value of conversion for particles of varying size is shown in Table 3.
|R||0.5 μm||1.0 μm||1.5 μm||2.0 μm||2.5 μm||3.0 μm|
|G||τ, s||τ, s||τ, s||τ, s||τ, s||τ, s|
Simulations R1 – R6 can be found in the Appendix of this work. The effect of particle radius on the required time to reach the desired conversion is shown in Figure 2.
Looking at the model (Equation 12) and Figure 2, it is obvious that a twofold increase in particle radius results in a fourfold increase in the required elapsed time to reach the same conversion, if all other conditions are unchanged.
The effect of K on the conversion is assessed in the additional 3 simulations (K1 – K3) – in each of these, K is set to a constant value – respectively, in the expected range the minimum, arithmetic mean and the maximum value is used here, based on the results of the S1 – S18 simulations. In all of these simulations the particle radius is fixed and is 1 μm. Simulations K1 – K3 can be found in the Appendix of this work.
The G = f(τ) relationship is shown in Figure 3, and is a clear representation of how nominal kinetic constant K influences the time required to reach a desired conversion, and is a graphic representation of which is also clear from the model Equation 12. Looking at the model (Equation 12) and Figure 3, it is clear that a twofold increase in nominal kinetic constant results in a twofold decrease in the required elapsed time to reach the same conversion, if all other conditions are unchanged.
In order to evaluate the mathematical model and to determine the effect of factors R and K on the conversion, the combined analysis is necessary, using descriptive statistics and dispersion analysis ANOVA, two-factor without replication .
The analyzed particle radius R and nominal kinetic coefficient K effect on conversion G = 0.90, using results of aforementioned simulations S1 – S18 and of the additional 18 simulations S19 – S36 in order to obtain a fuller picture of the influence. Results of these additional simulations can also be found in the Appendix.
Based on the results of the dispersion analysis shown in Table 6 it is deduced that the percentile effect of factor R is 95.38%, of factor K is only 2.80%, and the error scatter (unexplained dispersion) is 1.82%. However, based on the Fisher’s criterion at α = 0.05, the effect of both factors is significant, since FR = 105.040 > F0.05 = 3.326 and FK = 7.712 > F0.05 = 4.103. Since the effect of both factors R and K is significant and the numbers of residual degrees of freedom are νe, R = 5 < 10 and νe, K = 2 < 10, the dispersion analysis can be stopped here. Based on the analysis it is thus possible to conclude that the effect of the particle size on conversion is significantly larger than that of the nominal kinetic constant. These results indicate that the synthesis time during wet precipitation can be significantly decreased by decreasing the particle size.
3.1.2 Effect of conversion, particle size and time on the nominal kinetic constant value in the model
In the analysis a 2n factorial design, or Yates Analysis, is used. Complete factorial design (CFD) – a set of multiple measurements or simulations, which satisfies the following conditions: (1) The number of measurements or simulations is 2n, where n is the number of factors; (2) Each factor has only 2 values, the lowest and the highest; (3) During the measurement or simulation the lowest and the highest values are combined in all possible combinations. The advantage of such a design is the ease of parameter evaluation .
The influence of G, R and τ on the resulting value of K in the model is analyzed. 23=8 simulations are required in order to satisfy the aforementioned condition. The parameter configuration for the respective simulations is reported in Table 7, where „-” and „+” stands for the lowest and highest value, respectively.
Based on the earlier described value range the input parameter values are the following: Rmin = 0.5 μm, Gmin = 0.90, τmin = 18000 s, Rmax = 1.5 μm, Gmax = 0.95, τmax = 25200 s. Additional 8 simulations are performed (FD1 – FD8), based on the model equation in the form of Equation 16.
Results of the simulations can be found in the Appendix of this work. The results of the simulations are shown in the form of a matrix as represented in Table 7 according to the Yates analysis condition.
All 8 combinations and simulations are reported in the specific order as shown in Table 7. All responses are sorted as shown in Table 8. The values of the first column are then calculated – first 4 values are obtained by summing response pairs. For instance, the first value of the first column is the sum of a mean value and the effect’s response of R; the second value is obtained by summing G and R∙G responses. The last 4 values in the first column are calculated in an analogous fashion, but with the detraction of the responses. For example, the fifth value of the first column is found by detracting response mean from the effect’s response of R; the sixth value is found by detracting the G response value from the R∙G response.
|Gradation classes||K, m2/s, if τ = 5 h||K, m2/s, if τ = 6 h||K, m2/s, if τ = 7 h|
|R = 0.5 μm||2.563⋅10-18||2.136⋅10-18||1.831⋅10-18|
|R = 1.0 μm||1.025⋅10-17||8.544⋅10-18||7.324⋅10-18|
|R = 1.5 μm||2.307⋅10-17||1.923⋅10-17||1.648⋅10-17|
|R = 2.0 μm||4.101⋅10-17||3.418⋅10-17||2.930⋅10-17|
|R = 2.5 μm||6.408⋅10-17||5.340⋅10-17||4.577⋅10-17|
|R = 3.0 μm||9.228⋅10-17||7.690⋅10-17||6.591⋅10-17|
|R = 0,5 μm||3||6.53⋅10-18||2.18⋅10-18||1.35⋅10-37|
|R = 1,0 μm||3||2.61⋅10-17||8.71⋅10-18||2.17⋅10-36|
|R = 1,5 μm||3||5.88⋅10-17||1.96⋅10-17||1.10⋅10-35|
|R = 2,0 μm||3||1.04⋅10-16||3.48⋅10-17||3.46⋅10-35|
|R = 2,5 μm||3||1.63⋅10-16||5.44⋅10-17||8.46⋅10-35|
|R = 3,0 μm||3||2.35⋅10-16||7.84⋅10-17||1.75⋅10-34|
|K, m2/s, τ = 5 h||6||2.33⋅10-16||3.89⋅10-17||1.18⋅10-33|
|K, m2/s, τ = 6 h||6||1.94⋅10-16||3.24⋅10-17||8.17⋅10-34|
|K, m2/s, τ = 7 h||6||1.67⋅10-16||2.78⋅10-17||6.00⋅10-34|
|Scatter||Sum of Squares||Percentile effect||Degree of Freedom||Dispersion||Fisher’s criterion||p–value||F0.05,v,vZ|
|Effect of R||1.27⋅10-32||95.38||5||2.54⋅10-33||105.040||2.61⋅10-8||3.326|
|Effect of K||3.74⋅10-34||2.80||2||1.87⋅10-34||7.712||0.00941||4.103|
|Effect||Response||1st Column||2nd Column||3rd Column||Response Squares||Effect Value||Percentile Effect %||Effect’s (correlation) sign|
The second and third column is calculated in the similar manner as the first column, but in the place of responses first and second column values are used, respectively. This algorithm is then continued until the number of numbered columns is the same as the number of factors, which in this case is equal to 3. Values of the effect are then obtained by diving last column’s values with the number of simulation, meaning by 8 in this case. If the value is positive or negative, then the correlation sign is „+” and „-”, respectively. The values of effect are relative, but they allow to evaluate the percentile effect of parameters on the value of K. Percentile values are then obtained by dividing the modulus of effect values by the modulus of the sum of values .
Furthermore, it is possible to evaluate the effect sign; if a sign is „+” or „-”, then by increasing the parameter, the response increases or decreases, respectively. It can be concluded that the highest effect on the calculated value of the nominal kinetic constant K is the particle radius R, and is 59.95% as can be seen in Figure 4 and Table 8. The analysis also shows the coupled effect of parameters, which is most significant for the particle radius coupled with time. Thus, in order to obtain K values with a higher precision using the GB model, the highest attention has to be given to the measurements of the particle radius R.
In this short and simplistic preliminary study, a mathematical assessment of the Ginstling-Brounstein (GB) model showed that it can be potentially useful in describing the diffusion-controlled synthesis as in the case of wet precipitation (WP) of hydroxyapatite (HAp) using values of the parameters characteristic to the laboratory-scale synthesis of HAp. Particle size determination has the most significant effect (59.95%) on the kinetic coefficient value calculations and thus the highest attention has to be given to the careful measurements of the particle size in order to obtain K values with a high precision. It should also be noted, that reducing particle radius can significantly reduce diffusion-controlled reaction time in order to reach a desired conversion, even operating within the characteristic range of WP parameters in practice, and thus the overall synthesis time can be significantly reduced. An experimental validation of the model will follow.
The first author (Andrey) is thankful to doc. Ilo Dreijers (Ilo Dreyer) for continuous intellectual and moral support, as well as often highly exciting discussions on the topic and related subjects.
Conflict of interest: Authors state no conflict of interest.
 Qifeng S., Jiayun Z., Baijun Y., Jianhua L., Phase formation mechanism and kinetics in solid-state synthesis of undoped and calcium-doped lanthanum manganite, Mater. Res. Bull., 2009, 44(3), 649-653.10.1016/j.materresbull.2008.06.022Search in Google Scholar
 Tabero P., Bosacka M., Kurzawa M., The reaction mechanism and kinetics of the Zn2.5VMoO8 phase synthesis in the solid state, J. Therm. Anal. Calorim., 2001, 65(3), 865-869.10.1023/A:1011940502316Search in Google Scholar
 Boonchom B., Puttawong S., Thermodynamics and kinetics of the dehydration reaction of FePO4·2H2O, Physica B: Condens. Matter, 2010, 405(9), 2350-2355.10.1016/j.physb.2010.02.046Search in Google Scholar
 Khachani M., El Hamidi A., Kacimi M., Halim M., Arsalane S., Kinetic approach of multi-step thermal decomposition processes of iron(III) phosphate dihydrate FePO4•2H2O, Thermochim Acta, 2015, 610, 29-36.10.1016/j.tca.2015.04.020Search in Google Scholar
 Akbar J., Iqbal M.S., Massey S., Masih R., Kinetics and mechanism of thermal degradation of pentose- and hexose-based carbohydrate polymers, Carbohydr. Polym., 2012, 90(3), 1386-1393.10.1016/j.carbpol.2012.07.008Search in Google Scholar PubMed
 Gao P., Wang H., Jin Z., Study of oxidation properties and decomposition kinetics of three-dimensional (3-D) braided carbon fiber, Thermochim Acta, 2004, 414(1), 59-63.10.1016/j.tca.2003.11.017Search in Google Scholar
 Wionczyk B., Apostoluk W., Charewicz W.A., Adamski Z., Recovery of chromium(III) from wastes of uncolored chromium leathers. Part I. Kinetic studies on alkaline hydrolytic decomposition of the wastes, Sep. Purif. Technol., 2011, 81(2), 223-236.10.1016/j.seppur.2011.07.036Search in Google Scholar
 Ptáček P., Nosková M., Brandštetr J., Šoukal F., Opravil T., Mechanism and kinetics of wollastonite fibre dissolution in the aqueous solution of acetic acid, Powder Technol., 2011, 206(3), 338-344.10.1016/j.powtec.2010.09.040Search in Google Scholar
 Tudorachi N., Mustata F., Curing and thermal degradation of diglycidyl ether of bisphenol A epoxy resin crosslinked with natural hydroxy acids as environmentally friendly hardeners, Arabian J. Chem., (in press), 10.1016/j.arabjc.2017.07.008.Search in Google Scholar
 Corradini E., Pineda E.A.G., Hechenleitner A.A.W., Lignin-poly (vinyl alcohol) blends studied by thermal analysis, Polym. Degrad. Stab., 1999, 66(2), 199-208.10.1016/S0141-3910(99)00066-XSearch in Google Scholar
 Filimonov V.Yu., Koshelev K.B., Sytnikov A.A., Thermal modes of heterogeneous exothermic reactions. Solid-phase interaction, Combust. Flame, 2017, 185, 93-104.10.1016/j.combustflame.2017.06.020Search in Google Scholar
 Varganici C.-D., Marangoci N., Rosu L., Barbu-Mic C., Rosu D., Pinteala M., Simionescu B.C., TGA/DTA–FTIR–MS coupling as analytical tool for confirming inclusion complexes occurrence in supramolecular host–guest architectures, J. Anal. Appl. Pyrolysis, 2015, 115, 132-142.10.1016/j.jaap.2015.07.006Search in Google Scholar
 Blečić D., Ẑivković Ẑ.D., Martinović M.A, New method for the determination of reaction kinetics from DTA and TG curves. Part I. Definition of the method, Thermochim Acta, 1983, 60, 61-68.10.1016/0040-6031(93)85006-USearch in Google Scholar
 Ozola R., Krauklis A., Burlakovs J., Vincevica-Gaile Z., Rudovica V., Trubaca-Boginska A., Borovikova D., Bhatnagar A., Vircava I., Klavins M., Illite clay modified with hydroxyapatite– innovative perspectives for soil remediation from lead (II), Int. J.Agric. Environ. Res. (IJAER), 2017, 3(2), 177-189.Search in Google Scholar
 Sadat-Shojai M., Khorasani M.-T., Dinpanah-Khoshdargi E., Jamshidi A., Synthesis methods for nanosized hydroxyapatite with diverse structures, Acta Biomater., 2013, 9(8), 7591-7621.10.1016/j.actbio.2013.04.012Search in Google Scholar PubMed
 Dong L., Zhu Z., Qiu Y., Zhao J., Removal of lead from aqueous solution by hydroxyapatite/magnetite composite adsorbent, Chem. Eng. J., 2010, 165(3), 827-834.10.1016/j.cej.2010.10.027Search in Google Scholar
 Kanehira S., Kanamori S., Nagashima K., Saeki T., Visbal H., Fukui T., Hirao K., Controllable hydrogen release via aluminum powder corrosion in calcium hydroxide solutions, JAsCerS, 2013, 1(3), 296-303.10.1016/j.jascer.2013.08.001Search in Google Scholar
 Bērzina-Cimdina L., Dreijers I., Kreicbergs I., Hypothesis of Ca(OH)2 and H3PO4 reaction mechanism on solid particle surface, International Forum of Young Researchers “Topical Issues of Subsoil Usage”, St. Petersburg, Russia, 2011.Search in Google Scholar
 Grathwohl P., Diffusion in natural porous media: contaminant transport, sorption/desorption and dissolution kinetics, Kluwer Academic, Boston, 1998.10.1007/978-1-4615-5683-1Search in Google Scholar
 Fogler H.S., Elements of chemical reaction engineering, 4th ed., Pearson Education, Westford, 2005.Search in Google Scholar
 Renkin E.M., Filtration diffusion and molecular sieving through porous cellulose membranes, J. Gen. Physiol., 1954, 38(2), 225-243.Search in Google Scholar
 Thambynayagam R.K.M., The diffusion handbook: applied solutions for engineers, McGraw-Hill Professional, New York, 2011.Search in Google Scholar
 Ginstling A.M, Brounstein B.I., О диффузионной кинкTике реакций в сферических частицах, Журнал прикладной хмии, 1950, 23(12), 1249-1259, (In Russian).Search in Google Scholar
 Jander W., Reaktionen im festen Zustande bei höheren Temperaturen. Reaktionsgeschwindigkeiten endotherm verlaufender Umsetzungen, Zeitschr. für anorgan. und allgem. Chem., 1927, 163(1), 1-30, (In German).10.1002/zaac.19271630102Search in Google Scholar
 Arhipova I., Bālina S., Statistika ekonomikā. Risinājumi ar SPSS un Microsoft Excel. Mācību līdzeklis, 2nd ed., Datorzinību Centrs, Rīga, 2006, (In Latvian).Search in Google Scholar
 Box G.E.P., Hunter W.G., Hunter J.S., Statistics for experimenters: design, innovation, and discovery, 2nd ed., Wiley, Hoboken, 2005.Search in Google Scholar
© 2018 Andrey E.Krauklis, Ilo Dreyer, published by De Gruyter
This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License.