| Author |
 Topic  |
|
|
Karl-Ranz
Germany
Posts |
 Posted - 02/26/2008 : 1:27:36 PM
|
Origin 8.0 on Vista
Hello everyone
I am struggling for more than 2 weeks by now with a distribution fit and hope that someone in this forum might have experience with the stable symmetric Levy distribution. The problem I have is, that this distribution has no closed form, only an integral form.
I would like to fit the following distribution to a histogramm of data:
P(x,A,a,b)= A* 1/2pi * integral(-infinity..infinity) exp(-a*|k|^b) * cos(k*x) dk
where b is from 0 to 2 and A, a and b are the wanted parameters. This function has only two analytically solvable forms. For b = 2 we get a Normal-distribution and for b = 1 a Cauchy distribution.
As I have not yet found an comparable build-in function to use, I thought that creating this fit-function myself would be an good idea. Yet, I do not know how the handle the integral with the infinities and the abs() in the exp(). I am fairly new to the Origin C programming language and I do not know if routines for this problem exist. What I do know, is that the GSL C-Libraries have a routine for this distribution, but I am not sure if they can be used in Origin.
Online, you find many articles, where people have used a Levy fit, but nowhere is explained how they did it :) I think this fit is pretty common, yet a bit tricky to do.
Any help would be greatly appreciated.
Have a good evening and thank you very much.
Yours
Karl
Edited by - Karl-Ranz on 02/26/2008 1:39:09 PM
Edited by - Karl-Ranz on 02/26/2008 1:41:12 PM |
|
|
Echo_Chu
China
Posts |
Posted - 02/27/2008 : 04:03:57 AM
|
Please check this post to see whether it can help
Echo
OriginLab Corp |
 |
|
|
Karl-Ranz
Germany
Posts |
Posted - 02/27/2008 : 04:51:09 AM
|
Good morning
thank you very much for the quick reply. It looks like it describes the problem very well.
Thanks
Karl |
|
|
|
Karl-Ranz
Germany
Posts |
Posted - 02/27/2008 : 12:19:18 PM
|
Hello Echo_Chu and everyone else
The link you gave me was indeed helpful yet it only solved the problem, when we would integrate from 0...infinity, yet I unfortunately need the negative infinity as well and that is causing problems.
So far, based on the material, I have written this approach to get a hold on this problem:
-----------------------------------------------------------------
#include <math.h>
double dIntegral = 0.0; // integral value
double dInt = 0.0; // container for trapezium rule step
double dPrecision = 1e-6; // control integration to infinity
double k = 0.0; // 1st: Integration: 0..infinity
double dk = 0.05;
double f1, f2, g1, g2; // exp() and cos() for k and k+dk
do
{
f1 = exp( -SIGMA * k^ALPHA); // positive k => exp(-k)
f2 = exp( -SIGMA * (k+dk)^ALPHA);
g1 = cos( k*x );
g2 = cos( (k+dk)*x );
dInt = 0.5 * (f1*g1 + f2*g2) * dk; // Trapezium rule
dIntegral += dInt;
k += dk;
}while( fabs(dInt/dIntegral) > dPrecision ); // cos() could be negative => fabs(), header included
// DONE with 0..infinity, now 0..-infinity
k = -0.05; //reinitialize to go into negative direction
dk = -0.05;
do
{
f1 = exp( SIGMA * k^ALPHA); // k negative => exp(k), remember |k|
f2 = exp( SIGMA * (k+dk)^ALPHA);
g1 = cos( k*x );
g2 = cos( (k+dk)*x );
dInt = 0.5 * (f1*g1 + f2*g2) * dk;
dIntegral -= dInt; // "-=", as we are going from 0..-infinity
k += dk;
}while( fabs(dInt/dIntegral) > dPrecision );
y = Amplitude * dIntegral;
--------------------------------------------------------------
The parameters ALPHA, SIGMA, Amplitude get boundaries (Sorry, I did not stick with the parameter definition above: SIGMA==a and ALPHA==b):
0,5 < ALPHA < 2
0 < SIGMA < 100
0 < Amplitude < 1000
The initial parameter values are ALPHA=2, SIGMA=10 and Amplitude=500. The latter two relate fairly well to the data I see.
The problem with this program is, that the fit does not converge (almost as expected :) ) and Origin gets stuck when hitting (initialize parameters).
I would be glad if someone might have a look at the code to see what might be the matter.
The problem should also be of great interest to anyone in the finance sector, as the Levy distribution is used a lot in this field.
Thank you very much for any help.
Karl |
|
|
|
Echo_Chu
China
Posts |
Posted - 02/28/2008 : 05:01:06 AM
|
Hi, Karl
It looks like initial values of your function are not good. Would you mind to send your OPJ to tech@originlab.com and let's look into it?
Echo
OriginLab Corp |
|
|
|
Karl-Ranz
Germany
Posts |
Posted - 02/28/2008 : 6:35:38 PM
|
Hello Echo_Chu
Thank you for the offer, but unfortunately I cannot send this data.
I guess you are right in your assumption of wrong initial parameters. I have done runs in Excel, where I calculated the integral basically the same way as in Origin (using Visual Basic) and held it against the data. The best looking distributions were achieved with small values of SIGMA. Yet, the distribution was far away from describing the data as chi^2-tests showed.
Althought Excel gave an hint, a true converged fit would be better. In Excel I simply calculated the Levy Distribution for different values of ALPHA and SIGMA and tested each result with chi^2. Efficient, but far away from precise :)
So far, I have written a leaner code, which should do the same:
-----------------------------------------------------------------------------------------
double f1, f2, g1, g2;
double dInt = 0.0;
double dIntegral = 0.0;
double dPrecision = 1e-6;
double k = -10000.0;
double dk = 0.05;
do{
if(k < 0.0){
f1 = exp( SIGMA * k ^ ALPHA);
}
else f1 = exp( -SIGMA * k ^ ALPHA);
if( (k+dk) < 0.0 ){
f2 = exp( SIGMA * (k+dk) ^ ALPHA);
}
else f2 = exp( -SIGMA * (k+dk) ^ ALPHA);
g1 = cos( k * x );
g2 = cos ( (k+dk) * x );
dInt = 0.5 * (f1*g1 + f2*g2) * dk;
dIntegral += dInt;
k += dk;
}while( k < 10000.5 && fabs(dInt/dIntegral) > dPrecision );
y = Amplitude * dIntegral;
-----------------------------------------------------------------------------------------
What I would be interested in, is how Origin handles the parameters during a fit? Are they just simply looped within the boundaries, while the difference between fit and data is checked or are the handled differently? Knowing this could probably make the search for good initial parameters more easy. Sorry, for the inconvenience and thank you.
Karl
|
|
|
|
larry_lan
China
Posts |
Posted - 02/28/2008 : 10:10:17 PM
|
Hi Karl:
Origin use the Levenberg-Marquardt algorithm to fit nonlinear models. This link gives some outline about the algorithm.
Thanks
Larry
OriginLab Technical Services |
|
|
| |
Topic |
|
|
|
Stable Symmetric Levy Distribution Fit