1
$\begingroup$

I am using the code found here to fit multiple peaks of a large dataset which I don't know anything about. This was the only code I could find on the internet to do such a thing. So this is my attempt (it is pretty much the same code with the beginning deleted and a residual plot added) of a certain dataset and the results are good, it gets the correct number of peaks and it works.

However when trying it with a larger dataset, and the one I actually need to analyze, the code fails (it gives a lot of errors). I am guessing there is only a small technical error in the code or dataset somewhere. This is the dataset I need to analyse. That data was made by importing a text document and transposing it with integers.

The data:

data = Uncompress@Import["http://pastebin.com/raw.php?i=hkhuyvza"];

Peak function:

peakfunc[A_, μ_, σ_, x_] = A^2 E^(-((x - μ)^2/(2 σ^2)))

Model:

Clear[model]
model[data_, n_] := 
 Module[{dataconfig, modelfunc, objfunc, fitvar, fitres}, 
  dataconfig = {A[#], μ[#], σ[#]} & /@ Range[n];
  modelfunc = peakfunc[##, fitvar] & @@@ dataconfig // Total;
  objfunc = 
   Total[((Sqrt[data[[All, 2]]])/
       data[[All, 
         1]]) (data[[All, 2]] - (modelfunc /. fitvar -> # &) /@ 
         data[[All, 1]])^2];
  FindMinimum[objfunc, Flatten@dataconfig]]
$\endgroup$
  • $\begingroup$ We need to guess the errors you get? you need to show more effort making your question. Here its considered helpful and polite show you own efforts in the question and share your data and code attempts in a well formatted form, in the question, so we can quickly see the problem you are facing without going through the links and executing into Mathematica . Please help us to help you and edit your question accordingly. $\endgroup$ – rhermans Oct 30 '15 at 10:43
  • 1
    $\begingroup$ I have edited your question showing how to share your data and code properly. I did only part of the code, you do the rest. Please edit your question to improve it. $\endgroup$ – rhermans Oct 30 '15 at 10:48
2
$\begingroup$

Adapting @Silvia´s code (errr ... mostly copying it)

data = Uncompress@Import["http://pastebin.com/raw.php?i=hkhuyvza"];
data1 = Rest@Transpose[Rescale /@ (Transpose@data)];
peakfunc[A_, μ_, σ_, x_] = A^2 E^(-((x - μ)^2/(2 σ^2)));
Clear[model, modelvalue]
model[data_, n_] := Module[{dataconfig, modelfunc, objfunc, fitvar, fitres},
  dataconfig = {A[#], μ[#], σ[#]} & /@ Range[n];
  modelfunc = (peakfunc[##, fitvar] & @@@ dataconfig // Total);
  objfunc =  Total[((Sqrt[data[[All, 2]]])/ data[[All, 1]]) (data[[All, 2]] - (modelfunc /. fitvar -> # &) /@ 
         data[[All, 1]])^2];
  FindMinimum[objfunc, Join[{}, Flatten@dataconfig]]]
modelvalue[data_, n_] /; NumericQ[n] := If[n >= 1, model[data, n][[1]], 0]
fitres = ReleaseHold[ Hold[{Round[n], model[data1, Round[n]]}] /. 
    FindMinimum[modelvalue[data1, Round[n]], {n, 5},  Method -> "PrincipalAxis"][[2]]] // Quiet

resfunc = Join[peakfunc[A[#], μ[#], σ[#], x] & /@ Range[fitres[[1]]] /. fitres[[2, 2]], {}]

Show@{Plot[Evaluate[resfunc], {x, 0, 1}, 
   PlotStyle -> ({Directive[Dashed, Thick, 
         ColorData["Rainbow"][#]]} & /@  Rescale[Range[Length[resfunc]]]), 
         PlotRange -> All, Frame -> True, Axes -> False], 
   Plot[Evaluate[Total@resfunc], {x, 0, 1}, 
       PlotStyle -> Directive[Thick, Red], PlotRange -> All, 
       Frame -> True, Axes -> False], 
  Graphics[{PointSize[.003], Black, Point@data1}]}

Mathematica graphics

$\endgroup$
  • $\begingroup$ I have a question about this code. I am a newbie of Mathematica. Is it necessary to transpose and Rescale? $\endgroup$ – Clarine Jan 6 at 11:42

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.