# Create a Rayleigh distribution fitting experimental data

44 views (last 30 days)

Show older comments

DANIELE PASINI
on 16 Jan 2022

Answered: William Rose
on 18 Jan 2022

Hi, i'm trying to fit a rayleigh distribution to experimental data, but even if I've found the optimal parameter B for the distribution, it results in a completely different one.

I've tried using histfit (which works but I can't use in my assignment), makedist and the distributionFitter app.

Here's the code i'm using:

peaks=findpeaks(Data(:,i),'MinPeakHeight',0.01);

peaks=sort(peaks);

B=raylfit(peaks); %find parameter for Rayleigh distribution

ray=makedist('Rayleigh','B',B);

raydist=pdf(ray,peaks);

peaksdist=histcounts(peaks,1000,'Normalization','pdf'); %experimental distribution

hold on

x=linspace(0,max(peaks),length(raydist));

plot(x,raydist,'LineWidth',2)

When using the distributionFitter app I get this result instead:

But the parameter B used is the same in both (0.031).

How can I obtain the same distribution using the makedist function?

##### 2 Comments

### Accepted Answer

William Rose
on 17 Jan 2022

I plotted your data.

name='Cartel1';

data=load(name);

Data=data.data; %extract array 'data' from the structure 'data'

N=length(Data);

figure; subplot(311), plot(1:N,Data); xlabel('Sample'); title(name);

subplot(312), plot(1:5000,Data(1:5000)); xlabel('Sample');

subplot(313), plot(1:500,Data(1:500)); xlabel('Sample');

Now I understand why you want to model the distribution of peak heights. The data is a recording which shows bursts of activity. It is approximately symmetric about zero, like an electromyogram. Perhaps the peak height distribution will be useful for comparing recordings from different experiments.

Therefere let us find all the positive and negative peaks (i.e. peaks of -Data). Let us plot their histograms separately.

%find positive and negative peaks

peakpos=findpeaks(Data,'MinPeakHeight',0.01);

peakneg=findpeaks(-Data,'MinPeakHeight',0.01);

fprintf('Found %d positive peaks and %d negative peaks. \n',length(peakpos),length(peakneg));

%plot superimposed histograms of the positive and negative peak heights

figure;

h1p=histogram(peakpos,'FaceColor','g'); hold on

h1n=histogram(peakneg,'FaceColor','r');

legend('Positive peaks','Negative Peaks'); title('Histogram');

Since the histograms are very similar, we combine the sets into a single set of peak heights, :

%The histograms show that the positive and negative peaks appear to have

%the same distribution, so we combine them into a single set.

peaks=[peakpos;peakneg]; Npk=length(peaks);

Now fit the peak heights with a Rayleigh distribution.

bhat=raylfit(peaks); %sqrt of max.likelihood estimate of b^2

fprintf('Best fit Rayleigh parameter: bhat=%.4f.\n',bhat);

Plot the histogram of peak heights and the corresponding values expected from the best fit.

figure;

h2=histogram(peaks);

BinCenters=(h2.BinEdges(1:end-1)+h2.BinEdges(2:end))/2;

numExp=Npk*h2.BinWidth*raylpdf(BinCenters,bhat); %expected number at each bin

hold on; grid on

ylabel('Number of Occurences');

titlestr=sprintf('$Histogram, N=%d, \\hat{b}=%.4f$',N,bhat);

title(titlestr,'Interpreter','latex');

plot(BinCenters,numExp,'-r.'); hold off

Interesting. The peak height distribution is clearly not from a Rayleigh distribution.

The theoretical ratio of mean to standard deviation, for a random variable with a Rayleigh distribution, is

.

The ratio is 1.109 for peaks, which is not close to the theoretical value, if peaks has a Rayleigh distribution. This is more evidence for the non-Rayleighness of peaks.

##### 2 Comments

### More Answers (4)

the cyclist
on 16 Jan 2022

According to the documentation for raylfit, it takes the raw data as input. Is there a reason you are not simply using

B=raylfit(Data(:,i));

to find the fit parameter?

##### 5 Comments

William Rose
on 16 Jan 2022

Edited: William Rose
on 16 Jan 2022

@DANIELE PASINI, here is code with simulated data to demonstrate. As others have said, your data may not be appropriate for Rayleigh (whih is square root of the sum of two independent normal random variables). Your samples should be in vector x. Compute N=number of data points before computing numExp.

I did not use histfit, since you said it was not allowed. I did use raylfit(). If raylfit() is not allowed, then compute the maximum likelihood estimate for b^2, and take the square root:

bsqdhat=dot(x,x)/(2*N); . %MLE for b^2

The square toot of bsqdhat is not necessarily a MLE for b, and it may not be unbiased, but this seems to be what raylfit() does, because the result is the same to 4 decimal places, whenever I compare the estimates.

%RayleighFitTest.m

%WCR 20220116

%Generate simulated data to fit.

%Replace this section with your data.

b=0.031;

N=1000; %number of samples

x=b*sqrt(-2*log(rand(1,N))); %x=vector of samples whose distribution we wish to esitmate

%fit the data

bhat=raylfit(x); %Matlab's fit

bsqdhat=dot(x,x)/(2*N); %MLE for b^2:

fprintf('Best fit Rayleigh parameter: bhat=%.4f, sqrt(b^2hat)=%.4f\n',bhat,sqrt(bsqdhat));

%plot the histogram of the data, and the number expected from the fit

h=histogram(x);

BinCenters=(h.BinEdges(1:end-1)+h.BinEdges(2:end))/2;

numExp=N*h.BinWidth*raylpdf(BinCenters,b); %expected number at each bin

hold on; grid on

ylabel('Number of Occurences');

titlestr=sprintf('$Histogram, N=%d, \\hat{b}=%.4f$',N,bhat);

title(titlestr,'Interpreter','latex');

plot(BinCenters,numExp,'-bx');

Good luck!

##### 0 Comments

William Rose
on 18 Jan 2022

Edited: William Rose
on 18 Jan 2022

[I edited this comment to correct typographical errors.]

histogram(peaks) plots the number of occurences of values within the range of x-values for that bin. The expected number of counts for that bin is the integral of the probability density function over the bin, times the total number of samples.

Since each bin is reasonably narrow, we can approximate the integral with the midpoint value times the width:

where p(x)=probability density at x, and , and = bin width. Therefore

.

The command

raydist=pdf(ray,peaks);

does not give the correct result because it does not multiply by the number of samples, and does not multiply by the bin width, and because the argument peaks is not the correct argument to use. Instead of peaks, you should pass the x-values at which you want to evaluate the probability density function. And those desired x values are the bin center points. You can see and understand from my code how I computed the bin center points.

I hope your work goes well. What is the source of your signal? Seismogram? Biological recording?

William Rose
on 18 Jan 2022

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!