Good day to developers and users of Ifeffit, I'd like to get help with my question of statistical calculations in Ifeffit (I use the version embedded in Demeter, but in console mode or from my Fortran little program: think of Larch using but there are no some features and I'm not sure). The task is to try making Bayes statistics, in the sense that was published in Krappe et al. papers from 1999,2002,2004 years and the book of Bretthorst. I haven't found the software for it and my own tries to write it are not completed because of no experience. I saw something in papers about BFEFFIT, an extension for such statistics, but haven't found it in the web or any other information. May be somebody knows what stage is it on. I need it because of complicate model (4 spheres, may be splitting of spheres because it is many-phase system with dopant phase) and need of regularization in Tikhonov sense (parametra are correlated and I can not determine which constraint or model is better without statistics, and with nonconstrained parametra it simply will not fit) I understood that it can be emulated in Ifeffit with help of restraints. All thing is to: -minimize chi**2=(dat-model)**2/epsilon_k**2+A*sum(xi-x0i)**2, where A is regularizer, xi - parametra in fit end, x0i - a priori values of parametra. -calculate pair correlations and errors so that we can calculate the second derivatives and variance matrix of chi**2 determines in this sence. But I can not understand about how restraint goes in fit and how IFEFFIT calculates chi_square (not reduced) and errors. From manuals of FEFFIT and IFEFFIT I supposed that there goes: -LS minimisation of (dat-model)**2+restraint**2, may be normed, but it doesn't matter -chi_square is this value divided by eps_k**2, may be divided by points num -for errors found variance of chi_square defined above, but then error bars are scaled to chi_reduced to eliminate noise dependence -pair correlations are corr(i,j)/sqrt(corr(i,i)*corr(j,j) and doesn't depend of normalisations What is right? I need the above defined chi**2=(dat-model)**2/epsilon_k**2+A*sum(xi-x0i)**2 and its variance matrix, because it defines my statistics, optimal A,errors and 'projections' of parametra that I need to verify model. I use the fit in k space (not q) without Fourier transforms or k weight, because it is hard to do for this statistics. (Filtration of high fr3equency components I can do with rectangle windows if needed, I think, but it is not much bigger than noise) I thought (dat-model)**2 in k space is in real part for points from kmax to kmin for defined dat array(my_spectrum.chi for instance) and fit.chi array. But I calculated chi_square in hand from these arrays and there was no right number. for test of 1 sphere chi_square=8e-3, (dat-model)**2=3e-4, (dat-model)**2/noise**2=2,37e1, (dat-model)**2/noise**2=9.5e-2 and so on (regularizer part A(x-x0)**2=0.91). For instance. But if it calculates in real and imaginary part of Fourier transform no meaning what space I use - and this function is used to find variance - then I'm not sure how to find error bars and variance matrix I need. There are many other questions and limitations, I want to make this calculation approximately with what I can (constaint noise, may be +bkg part as Ifeffit gives, two or more spectra - and the other things later), but I hope someone could help the beginner to solve this problem or find the simplest way in existing software... I need to know what is the stage of this topic and need I end my program for it or not. There in mind are such things as mu0 polynome fit, many-spectra fit with PDF and diffraction spectra for model and so on, but may be there is the best program in one step from completeness and I need not to do it. Thanks sincerely, Olga Kashurnikova (from MePHI XAFS group in Moscow)
Hi Olga,
It's a bit hard for me to tell what here is actually a question and
what is an assertion of preconceived ideas, or a description of what
you want (or think you want). You might find it more useful to ask
questions in smaller units, and aim for more clarity.
On Wed, Mar 4, 2015 at 3:51 PM, Ольга Кашурникова
Good day to developers and users of Ifeffit, I'd like to get help with my question of statistical calculations in Ifeffit
(I use the version embedded in Demeter, but in console mode or from my Fortran little program: think of Larch using but there are no some features and I'm not sure).
Not sure what that means. FWIW, there is not going to be future development or changes beyond critical bugs to Ifeffit. If you want to use anything supported or ready for adding features, Larch is the only way to go.
The task is to try making Bayes statistics, in the sense that was published in Krappe et al. papers from 1999,2002,2004 years and the book of Bretthorst.
Not familiar with that book.
I haven't found the software for it and my own tries to write it are not completed because of no experience.
I believe that Krappe and Rossner's used a software modified from Feffit, but I'm not sure they ever made this available for others to use.
I saw something in papers about BFEFFIT, an extension for such statistics, but haven't found it in the web or any other information. May be somebody knows what stage is it on.
I've never heard of this.
I need it because of complicate model (4 spheres, may be splitting of spheres because it is many-phase system with dopant phase) and need of regularization in Tikhonov sense (parametra are correlated and I can not determine which constraint or model is better without statistics, and with nonconstrained parametra it simply will not fit)
Tikhnonov regularization and Bayesian analysis seem quite different to me. Can you explain what you mean? Of course parameters are correlated... this happens all the time, but does not mean they cannot be fitted. Also, saying "simply will not fit" gives no clue about what might be happening.
I understood that it can be emulated in Ifeffit with help of restraints. All thing is to: -minimize chi**2=(dat-model)**2/epsilon_k**2+A*sum(xi-x0i)**2, where A is regularizer, xi - parametra in fit end, x0i - a priori values of parametra. -calculate pair correlations and errors so that we can calculate the second derivatives and variance matrix of chi**2 determines in this sence.
Restraints are a way to encode prior knowlege, but I'm not sure it's really a Bayesian analysis. Or perhaps: it's one way of doing a Bayesian analysis, but almost certainly not the only way. For a Bayesian analysis you need to find the most likely values for parameters p given data (and a model) and a set of prior probabilities for the parameters. So you need to: a) have an expectation value and probability distribution for each parameter b) somehow encode this in the actual analysis. Note that this doesn't say anything about how one does the analysis. When doing a Bayesian analysis, one is given an additional problem (identify a probability distribution for each parameter) and no guidance of what to do once that is known. In this sense, Bayesian analysis is "obviously better" because it prescribes a general maxim ("include prior knowledge of your parameters") and no guidance of any practical utility. It is approximately equivalent to saying "Find a career that makes you happy". Chi-square minimization is completely different -- a separate topic really. It's not Bayesian or Non-Bayesian (Frequentist), it's just a method. Practical, no theoretical. It finds a set of parameters p to make "[data - model(p)]/epsilon" match as well as possible. This is convenient because there are a variety of mathematical methods ("Curve Fitting / Optimization") for finding a set of parameters that minimizes some array in the least-squares sense. If looking at it from a theoretical viewpoint, this is exactly asserting that the variance in the data is given by a Normal / Gaussian distribution. If you don't know what probability distribution to use, this would be one of the first to assume. Given this, one route to a practical way to encode "prior information" / Bayesian priors is to restraints. Here one simply makes up a "data point" and extends the array to be minimized: [data-model(p)]/epsilon + restraint1 + restraint2 + .... If one sets restraint1 to be [paramA - A_0]/eps_A that would effectively assert a normal distribution for parameter A centered at A_0 with sigma = eps_A. Of course, that doesn't address the question of how much you believe "the parameters should be varied so that the model matches the data" compared to "the value for parameter A should be near A_0". You may need to adjust a weighting factor, or scale eps_A accordingly. That's sort of like regularization.... To be clear, Bayes gives no guidance on how to do this.
But I can not understand about how restraint goes in fit and how IFEFFIT calculates chi_square (not reduced) and errors. From manuals of FEFFIT and IFEFFIT I supposed that there goes: -LS minimisation of (dat-model)**2+restraint**2, may be normed, but it doesn't matter -chi_square is this value divided by eps_k**2, may be divided by points num -for errors found variance of chi_square defined above, but then error bars are scaled to chi_reduced to eliminate noise dependence -pair correlations are corr(i,j)/sqrt(corr(i,i)*corr(j,j) and doesn't depend of normalisations
See above. Chi-square is the scaled sum of squares of the array to be minimized. So, yes, it would include restraints. And restraints can alter the uncertainties and the correlations because they change the array to be minimized.
What is right? I need the above defined chi**2=(dat-model)**2/epsilon_k**2+A*sum(xi-x0i)**2 and its variance matrix, because it defines my statistics, optimal A,errors and 'projections' of parametra that I need to verify model.
I use the fit in k space (not q) without Fourier transforms or k weight, because it is hard to do for this statistics.
Sorry, I can't read this in any way except "I do X because it is hard to do". Is that really what you mean? I strongly recommend not fitting is k-space. Doing so willfully ignores something that you do, in fact, know about the model you are using for the data -- that it is limited in R-space. That is, if the model has no weight above R>8 Angstroms, why try to fit this part of the data? There is no situation in which fitting in k-space is better than fitting in R-space, and very many situations where it is worse.
(Filtration of high fr3equency components I can do with rectangle windows if needed, I think, but it is not much bigger than noise) I thought (dat-model)**2 in k space is in real part for points from kmax to kmin for defined dat array(my_spectrum.chi for instance) and fit.chi array. But I calculated chi_square in hand from these arrays and there was no right number. for test of 1 sphere chi_square=8e-3, (dat-model)**2=3e-4, (dat-model)**2/noise**2=2,37e1, (dat-model)**2/noise**2=9.5e-2 and so on (regularizer part A(x-x0)**2=0.91). For instance.
Sorry, I'm not sure what you are trying to say. Is that your by-hand calculation of chi-square is not matching the one reported? If you wanted us to help you solve this, you would have to give us more details.
But if it calculates in real and imaginary part of Fourier transform no meaning what space I use - and this function is used to find variance - then I'm not sure how to find error bars and variance matrix I need.
Not sure what you're talking about here.
There are many other questions and limitations, I want to make this calculation approximately with what I can (constaint noise, may be +bkg part as Ifeffit gives, two or more spectra - and the other things later), but I hope someone could help the beginner to solve this problem or find the simplest way in existing software... I need to know what is the stage of this topic and need I end my program for it or not. There in mind are such things as mu0 polynome fit, many-spectra fit with PDF and diffraction spectra for model and so on, but may be there is the best program in one step from completeness and I need not to do it.
Sorry, I don't understand this.
Thanks sincerely, Olga Kashurnikova (from MePHI XAFS group in Moscow)
I'm not sure I answered any of your questions, but hope it helps. Cheers, --Matt
Hello, Dr Matt Newville, It is the second part about Bayesian analysis itself. I hope I answer in a thread, not only to you. May be it isn’t in a thread, I don’t know how to mend it now The question was splitted because I was not sure if there are people who already made this thing with Bayesian statistics in program and I should ask about is it right at all before doing. The same is for the last paragraph of message you answered, I tried to say that there are some things I can’t do that simple with IFEFFIT and maybe they are worth of work even if the simple thing is already done. I was not sure I should make multiple threads for it because it is interconnected with the main question. For this part your answer is very informative and helpful. I thought of using Larch, but now there are not many new features for my task, and I didn’t find some old, such as mu0 co-fitting with XAFS function. Ifeffit does it for a fixed num of points, I think I need to vary it, but don’t know if it is possible. I will surely try to study how to use Larch but have seen it not a long ago and can’t say much of it now. I am thankful to you for the info it will be more extensible. I’d like to ask about it when I will learn more,if possible. About what is Bayes for. I have the complicated system with nanodomains of fluorite and pyrochlore phase. The samples were grown in special conditions and the metastable phases are practically valuable in this case, that’s why we should understand how the structure evolves. The task was to find the structure from XAFS spectra, but there should be the complicated model for exact characterization. Diffraction gave the evidence of multicomponent system with evolving proportion (pyrochlore domains grow inside fluorite). One of the structure features is the splitting of oxygen coordination sphere, that have to be in pyrochlore and not to be in fluorite. In the same time, fluorite structure is defected in the sense of dopant structure – two cation elements are mixed in fcc sublattice and cation-oxygen distance is different for them, so we can get the splitting of oxygen sphere to 5 spheres, for instance, but they will be undistinguishable from the wide Gaussian sphere I think, only in cluster model with constraints. The pyrochlore splitting is more evident, but I think it needs statistical proving. Also if I could make the more complicated multicomponent-with-dopant model (I tried to do it, but in the different sense, because we didn’t have the full anomalous diffraction data that time), it should be proved to give the better results. Some of this information is in correlations, but they can’t give the information of what parametra exactly should be eliminated from the model. And projections can. The book of G.Larry Bretthorst ‘Bayesian Spectrum Analysis and Parameter Estimation’ from the series ‘Lecture notes in Statistics’ (Springer-Verlag,1988) is not very different by concept from Krappe and Rossner papers which you seem to be familiar to. It is not about XAFS but about different data modeling. Author makes many tests with sum of one or many frequency sine functions, with or without broadening, with trend lines and so on, trying to calculate it more theoretically, but with a little program that can be used only after minimization. He proves the thing that two frequencies can be statistically distinguished, which model is better, with one or two frequencies. I think it is very close to sphere splitting. Also there is a formula for relative probability of two models which seems to be useful (he proves it can distinguish Gauss and Lorentz broadening), and there is a method to define which number of polynomes should be in the trend line (as mu0 – and may be it works for knots num too, I tried to check it, but it needs co-fitting of mu0 and chi). I thought something like this can help say the more definite words about the structure of my samples, that’s why I tried it. For XAFS spectra, I have chosen the Krappe and Rossner method of calculus now because it is difficult to use tabulated amplitudes and phases in the different method, but they are interconnected. What of the k-space, the math in Bayesian method is for minimization of (dat-model)^2/eps^2 itself, not FT of it. And statistical values (covariance matrix, projections et all) are calculated for it. Bretthorst says that FT works only for clear sine functions and gives the same as his analysis. But for complicated cases he insists on other analysis. But of course there is need to filter the high-frequency part from XAFS spectra, but I am not sure how to make it in Bayesian sense, and how to calculate errors, correlations and projections, and model comparison for this case. I thought I could filter high-frequency with rectangle windows (possibly with *k^2 and the result /k^2 because of XAFS formula) and will not lose any signal but the non-necessary high-frequency. And then to model it as usual. But of course it could be more right to use part of FT and to transform math, but I am not sure how it should be done now. I wanted to fit not-filtered and filtered functions and compare, because the high-frequency part is not much more than noise level and may it could be neglected in our case. I think the model takes into account that we have limited the R spectrum – the model consists only of the first coordination spheres? The k grid for this method needs not to be even, and the step is taken into account by the math of Bayesian method, or I mistake? The N_idp information in this case, I think, is given by parametra projections and other definition of “what part of parametra space is defined by spectra and what of a priori information”©Krappe and Rossner, I think? I said ‘is hard to do’ Bayesian analysis with Fourier-filtered spectrum (may be missed a word) because I didn’t find in literature how it should be done, and couldn’t solve this problem myself. I’m not sure I have the experience for it, though I thought of it and will search the solution. The method of statistics IFEFFIT use and Bayesian statistics are somewhat different in concept I thought, and I can’t use it that simply here, but I can mistake. May be there are books I didn’t find on this stage, may be the chi_square for transformed signal can work as well? I’m a novice (that has found a useful algorythm in papers) and am not sure in it, but may be it can’t be done all at once, that’s why I try to simplify the task to get some valuable structure information and apply an existing method first. May be it will answer why I need k-space and all these renormalizations and extra work. Thankful for your assistance, Olga Kashurnikova, MEPhI, Moscow
Hello, There are text files in attachment, save them to the same directory an load check.iff in Ifeffit shell I hope attachments will work well. There are two .iff files with macros I need and the retranslation of what program does. And feff and spectrum files, too. The spectrum is generated from the same feff, for 1 shell, but in VIPER for easy visualization, so it can be the other formula, but should be close (there were close parametra e1=5.1523,n1=7.0612,r1=2.1107,s1=7.4408e-3 - it is the fitted model for mu0 corrections initially). I haven't append the noise, think that foe chi_square formula check it is not important, but if it is, I can try. May be some other parametra wrong, I can't understand it now. Compare program chi_square and chi_square_calc that I calculated for the suggestion of formula. In this case the restraint=0, so the restraint part doesn't mean - it should be checked when we understand the normalization. Here, for 1 sphere I need not a regularizer, but the method can be tested on it and should also give Bayesian analysis good, because there are some correlations. (I tried it in the own program - it works, but there are other errors why I need to use IFEFFIT while haven't completed it) So, this formula doesn't work, what is wrong - how to take FT into account or renormalize? If I know how it renormalizes, I can renormalize the errors for my need. I understand that errors and correlations are of a covariance matrix calculated for chi_square(param)/chi_reduced(min), where the denominator is the end value in min and only a norm. I need a covariance of chi_square without N_idp/N_data multiplier Thankful, Olga Kashurnikova, MEPhI, Moscow
participants (3)
-
Matt Newville -
Olga Kashurnikova -
Ольга Кашурникова