Posts Tagged ‘Mathematica’

Mathematica in Examples – Replace all Elements in an Array that satisfy a Condition

Here is a short but useful code snippet, that replaces all elements in an array that satisfy a condition. Mathematica has many ways to do such things, but provides remarkable shortcuts.

Specifically,  I want all elements in the array

array={0.601, 0.9, 1.553, 2.588}

that are larger 1 to be set to 1 (thus clipping of all elements larger 1).

Download Mathematica Notebook: ReplaceElementsinanArray.nb

The fastest way to do that is to use the ReplaceAll function combined with the Condition function. First is highlighted in blue and latter in red.

In:=

array/.x_/;>1->1

The ReplaceAll function is given by /.condition->replacement , where x_/;>1 is our condition number greater 1 and the replacement “1″. The output gives than

Out:=

{0.601, 0.9, 1, 1}

Pretty simple, huh?

Alternatively, the same result is achieved by using the Table function together with the If function:

In:=

Table[If[array[[i]] > 1, 1, array[[i]]], {i, 1, Length[array]}]

The code is not as compact as the one above, but has the same output.

Mathematica in Examples – Fit Lorentz distribution with Constraints

Mathematica’s NonLinearModelFit function is useful to fit an arbitrary function to a set of data. In this example code a Lorentz function with constraints for the parameters is fit to some example data. Initial fit parameters are guessed.

File: nonlinearfit_lorentz_osfight.nb

1. We define the Lorentz function with linear background in Mathematica as

Lor[{A_,μ_,Γ,C0_,m_},ν_]:= (A Γ)/((ν-μ)^2+(Γ/2)^2) + m ν + C0;

where 4A/Γ is the peak height, μ the center, Γ the linewidth of the of the Lorentzian (the code can be copied to Mathematica). mν + C0 is the linear background and can be set to zero if not needed.

2. Lets take some random data and try to guess the initial parameters by plotting the guess function and the data together. The guessed initial values are colored in red

Input=

data = {{0.`, 0.0007568602316023421`}, {0.4`,
   0.0006077641663579275`}, {0.8`,
   0.0010957089335075778`}, {1.2000000000000002`,
   0.0014422342777222481`}, {1.6`, 0.0021026269660642818`}, {2.`,
   0.002035918694285509`}, {2.4000000000000004`,
   0.0019560986865232773`}, {2.8000000000000003`,
   0.002807918106577903`}, {3.2`, 0.0032891006273276686`}, {3.6`,
   0.009439170796893468`}, {4.`, 0.012570566138374808`}, {4.4`,
   0.01820541122261653`}, {4.800000000000001`,
   0.01555451311477995`}, {5.2`,
   0.022197008830661266`}, {5.6000000000000005`,
   0.013782511728324743`}, {6.`, 0.013581988241044292`}, {6.4`,
   0.0090081109450604`}, {6.800000000000001`,
   0.003354272775540456`}, {7.2`,
   0.0021712171237978674`}, {7.6000000000000005`,
   0.0030934291572321778`}, {8.`, 0.0019866310063820377`}, {8.4`,
   0.0010959725673603206`}, {8.8`,
   0.0012691490482224404`}, {9.200000000000001`,
   0.0006927858593134837`}, {9.600000000000001`,
   0.0011999520986368984`}, {10.`, 0.0007696824792981429`}};

guessPlot = Plot[Lor[{0.005, 5, 1, 0, 0}, x], {x, 0, 10}, {PlotStyle -> Red}, PlotRange -> All];
dataPlot = ListPlot[data];
Show[{guessPlot, dataPlot}, PlotRange -> All]

Output=

The guessed Lorentz curve (red) is shown next to the data points.

3. Our inital values look decent, so we can try to fit a curve to the data. The initial values from above are labeled in red and constraints in blue.

Input=

nlmLor=NonlinearModelFit[data,{Lor[{A,μ,Γ,m,C0},x],{Γ>0,A>0}},{{A,0.005},{μ,5},{Γ,1},{m,0},{C0,0}},x]; 
nlmPar=nlmLor["BestFitParameters"]

Output=

{A->0.00987806,μ->5.05304,Γ->2.18077,m->0.000195088,C0->-0.0000407613}

4. Now we can fit the result and check the fit quality
Input=

fitPlot = Plot[Lor[{A, μ, Γ, m, C0}, x] /. nlmPar, {x, 0,10}, PlotStyle -> {Blue, Thick}, PlotRange -> All];
nonlinearfit_lorentz_osfightShow[dataPlot, fitPlot]

Output=

The fit Lorentz curve (blue) is shown with the data points.

Tips

  • the “quality” of the initial parameters are quite important for the convergence of the fit
  • using constraints requires more computation effort and is not recommended with the NonLinearModelFit
  • detailed fit errors are given using for example nlmLor["ParameterTable"] or nlmLor["ParameterTableEntries"]

This article is part of a series explaining how to use Mathematica for experimental data analysis. In this series:

Part 1 – Fit Lorentz distribution with Constraints

Part 2 – Fit a triplet Gaussian (coming soon)

Return top