3
$\begingroup$

I would like to know how to plot a make a comparison between a histogram and a set of random data, I have the following probability density function (PDF) which is given for this problem. I can generate the histogram but I do not know how to add the PDF function and make the plot both together.

I have seen some question and problem using RandomVariate function but it does not work. This is the code.

Clear["Global`*"];  Remove["Global`*"];
 L = 10; Ns = 10^4;
 x = Table[RandomReal[L], {x, 0,Ns}];  
 y = Table[RandomReal[L], {y, 0, Ns}];
 z = Table[RandomReal[L], {z, 0, Ns}]; 
 sf = Abs[Sin[x^2  + y^2] + Cos[x y z]];
 Histogram[sf, {0.25}, "Probability"]

The (PDF) is f = k^2/(k^2 + 1)^(5/2) The goal is to obtain a similar plot like this, enter image description here , where the line represent in this case f.

$\endgroup$
  • $\begingroup$ There are a number of problems with your formulation. First, you do not need x = Table[RandomReal[L], {x, 0,Ns}]; the x on the RHS is superfluous, and likely to cause problems. All you need is: x = RandomReal[10, 10^4]; Second, your random data sf is constrained to lie between 0 and 2 (i.e. is bounded). Your pdf is not bounded (or if it is, it is not stated). Third, even if your density was unbounded, it does not integrate to unity, so is not well-defined. All that leaves aside the fact that the density does not fit the data. $\endgroup$ – wolfies Aug 19 at 15:15
  • $\begingroup$ @wolfies, Your are right all of your points, in the problem the PDF is non-normalizable, and the data is not bounded. I would really like to know how can I get the right PDF to approximate this data? Can you please provide an example. $\endgroup$ – irondonio Aug 23 at 9:21
  • $\begingroup$ I can't tell from your response if you're agreeing or disagreeing with @wolfies comment about the distribution being bounded: it is bounded by 0 and 2. I'm not convinced that there would be a simple form for the density but you could certainly estimate it with a nonparametric density estimate: skd = SmoothKernelDistribution[sf, Automatic, {"Bounded", {0, 2}, "Gaussian"}]; Show[Histogram[sf, {0.05}, "PDF"], Plot[{(15 Sqrt[5]/8) k^2/(k^2 + 1)^(5/2), PDF[skd, k]}, {k, 0, 2}]] . $\endgroup$ – JimB Aug 23 at 13:07
  • $\begingroup$ @JimB, I used your code and it works. I wonder if using SmoothKernelDistribution, What happens if we do not use "Gaussian" or another. I plot without this option and it also works. (For this case which function is used) It is important to know what is the analytical form the function ( plotted in blue color) given which is approximated by the Histogram? It is possible to obtain this function from Mathematica. or It is possible to build this function (PDF) from the data? $\endgroup$ – irondonio Aug 23 at 17:22
  • $\begingroup$ From a scientific point of view I'd recommend comparing CDFs, not PDFs. $\endgroup$ – Roman Aug 24 at 17:34
5
$\begingroup$

Use Show:

Show[
Histogram[sf, {0.25}, "Probability"],
Plot[f, {k, 0, 2}]
]

enter image description here

$\endgroup$
  • $\begingroup$ I think you want "PDF" rather than "Probability". $\endgroup$ – JimB Aug 19 at 13:53
  • $\begingroup$ You're right, I coppied/pasted the code from the OP. $\endgroup$ – Themis Aug 19 at 16:40
  • $\begingroup$ @JimB, Your comment is correct. I have a question If I have some histogram what is the way to obtain the right PDF to approximate this data? If you no have Idea about the data follow a particular distribution, what is the procedure and step to follow in Mathematica? $\endgroup$ – irondonio Aug 23 at 9:29
  • 1
    $\begingroup$ @irondonio See my comment on your original question above. Nonparametric density estimators should be used in place of histograms in almost all situations. $\endgroup$ – JimB Aug 23 at 13:14
2
$\begingroup$

The OP writes: "It is important to know what is the analytical form of the function (plotted in blue color) given, which is approximated by the Histogram"

It might be difficult to find the analytical form of sf. But we can approximate the analytical form quite easily by noting that your data appears exponentially downward sloping, but bounded on (0,2). So a natural choice would be a truncated $\text{Gamma}(a,b)$ model, truncated above at 2. This model will have pdf

$$g(s) = c s^{a-1} \exp{(-\frac{s}{b})}$$

where constant c = 1/(b^a*Gamma[a]) - Gamma[a, 2/b].

Parameters $a$ and $b$ can then be found using maximum likelihood.

Here is a plot of your histogram and the fitted truncated Gamma model when $a = 0.805$ and $b = 4.47$: enter image description here

...which seems to work nicely. The fitted density is then approximately something simple like:

0.56 Exp[-0.224 s] s^(-.2)
$\endgroup$
0
$\begingroup$

This is a specific function which can be modified for any PDF you like. (It looks horrible because i wrote it.) Set PrintDistributionAndFitOpt -> True to see the plot. Otherwise it just returns the fit.

Clear[GetGaussianFitToDataDistribution];
Options[GetGaussianFitToDataDistribution] = {
UseNumberBinsOpt -> 50
    ,PrintDistributionAndFitOpt -> False
};
GetGaussianFitToDataDistribution[InData_List, OptionsPattern[]] := Module[{
    UseNumberBins, PrintDistributionAndFit
    , DataScore, MaxBin, MinBin, BinWidth, UseBins
    , wsigma, wmu, x, wFit
},
{UseNumberBins, PrintDistributionAndFit} = OptionValue[{UseNumberBinsOpt, PrintDistributionAndFitOpt}];
If[EvenQ[UseNumberBins], UseNumberBins++;];
If[Length[InData] < UseNumberBins
, Print["Warning in GetGaussianFitToDataDistribution: Less data than bins."];
];
(* Determine the bins *)
MaxBin = Ceiling[Max[InData]];
MinBin = Floor[Min[InData]];
BinWidth = N[(MaxBin - MinBin)/UseNumberBins];
UseBins = Range[MinBin + 0.5*BinWidth, MaxBin, BinWidth];
(* Generate data to be fit using BinCounts
- normalize to 1.0 or the optimization gets confused
*)
DataScore = Transpose[{UseBins
, (N[#] &) /@ BinCounts[InData, {MinBin, MaxBin, BinWidth}]
}];
DataScore = ({#[[1]], #[[2]]/Max[DataScore[[All, 2]]]} &) /@ 
DataScore;
(*Print[DataScore];*)
(* Find the Best Fit 
- Uses Mean and StDev for estimates of the partameters
*)
wFit = NonlinearModelFit[DataScore
, PDF[NormalDistribution[wmu, wsigma], x]/
PDF[NormalDistribution[wmu, wsigma], wmu]
, {{wmu, Mean[InData]}
, {wsigma, StandardDeviation[InData]}}, x];
If[PrintDistributionAndFit, Print[Show[{
ListPlot[DataScore, PlotStyle -> Darker[Red], PlotRange -> All
, PlotLegends -> Placed[{"Original Data"}, {0.85, 0.85}]
]
, Plot[wFit[x], {x, MinBin, MaxBin}
, PlotLegends -> Placed[{"Fitted Gaussian"}, {0.85, 0.75}]
]
}
, PlotLabel -> ToString[{wmu, wsigma} /. wFit["BestFitParameters"]]
    <> "\t Samples = " <> ToString[Length[InData]]
, FrameLabel -> {"Bin Mid-Line", "Relative Frequency"}
]];
];
{wmu, wsigma} /. wFit["BestFitParameters"]
];
(*
TestData=RandomVariate[NormalDistribution[1000.,20.],10^5];
GetGaussianFitToDataDistribution[TestData
,PrintDistributionAndFitOpt\[Rule]True
,UseNumberBinsOpt\[Rule]50
]
*)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.