Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to set up a system of 400*N equations to deconvolute spectral data for four different compounds.

I monitored 400 different wavelengths in N solvents (3 solvents in this case). In order to do this, I want to try to minimize the following function:

$A_{\lambda_i, N_i} - (\epsilon_{1 (\lambda_i, N_i)}C_1 + \epsilon_{2 (\lambda_i, N_i)}C_2 + \epsilon_{3 (\lambda_i, N_i)}C_3 + \epsilon_{4 (\lambda_i, N_i)}C_4 = 0 $

Here $A$ is the absorbance of the spectrum I want to deconvolute at specified wavelength ($\lambda$) in solvent ($N$).
$\epsilon$ is the molar extinction coefficient for each different compound that corresponds to the absorbance at the given wavelength and solvent. Then unknown quantities are $C_1$ to $C_4$

https://drive.google.com/folderview?id=0B5gGZ4-SD7-ZdHFveVY5b2J3OTQ&usp=sharing

My data is organized as follows.

Each row has 5 entries and I have 400 rows one (for each wavelength) [$A \; \epsilon_1 \; \epsilon_2 \; \epsilon_3 \; \epsilon_4$]

Each solvent has a table organized just like this one. How do I used mathematica to set up this system of equations and get the error associated with the minimization? Also, the solution needs to be constrained such that C$_n$ is positive.

Import the data (uploaded on google drive)

phosBB = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_BB.CSV"];
phosMR = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_MR.CSV"];
phosPT = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_PT.CSV"];
phosTB = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_TB.CSV"];
phosUI = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_UI.CSV"];

PreProcess data to get rid of headings:

phosBB = Delete[phosBB, 1];
phosMR = Delete[phosMR, 1];
phosPT = Delete[phosPT, 1];
phosTB = Delete[phosTB, 1];
phosUI = Delete[phosUI, 1];

Calculate $\epsilon$ from the known solutions (BB-TB):

\[Epsilon]PhosBB = phosBB[[All, 2]]/concPhosBB;
\[Epsilon]PhosMR = phosMR[[All, 2]]/concPhosMR;
\[Epsilon]PhosPT = phosPT[[All, 2]]/concPhosPT;
\[Epsilon]PhosTB = phosTB[[All, 2]]/concPhosTB;

Setup data table [$A_{unk} \; \epsilon_1 \; \epsilon_2 \; \epsilon_3 \; \epsilon_4$]

phosDataTable = 
  Transpose[{phosUI[[All, 1]], \[Epsilon]PhosBB, \[Epsilon]PhosMR, \[Epsilon]PhosPT, \[Epsilon]PhosTB}];
share|improve this question
1  
Have you attempted something already with e.g. FindFit? Can you show your code? –  Mike Honeychurch Feb 22 at 2:37
3  
It looks like a linear system -- so you can rewrite your equations in matrix form m.x==b and then use LinearSolve. m is your data {A, e1, e2, e3, e4}, b is zero and x={c1, c2, c3, c4}. –  bill s Feb 22 at 3:43
    
@MikeHoneychurch I'm not sure what mathematica function to use. The idea is to guess C1-C4 to try to get as close to zero as possible. There will be some error associated with the calculation because I don't think there is any exact solution. we just want to get close. –  olliepower Feb 22 at 21:14
    
@bills I tried using Linear solve where m was {e1,e2,e3,e4}, x = {c1,c2,c3,c4}, b = {A$_{unk}}. (a rearrangement of the equation i posted above). There was no exact solution. –  olliepower Feb 22 at 21:16
2  
400 points is overkill for visible spectra of four components. If you have access to Harris' Quantitative Chemical Analysis, chapter 18 provides a very thorough explanation and demonstrates the solution using Excel (which turns out to be a MMA one or two liner). –  bobthechemist Feb 23 at 0:21

2 Answers 2

up vote 4 down vote accepted

Import the data and filenames, cull the datasets of information that is not useful:

data = Import /@ FileNames["*.csv"];
filenames = FileNames["*.csv"];
data = data[[All, 2 ;;, {1, 2}]];
ListLinePlot[data, PlotLegends -> filenames]

Mathematica graphics

Looks like PT data was either incorrectly collected or is effectively transparent in the region of interest. The question currently is missing concentration information, so I will assume each standard is 0.1 M just for pedagogical purposes. Convert the standard ordinates to molar absorptivities.

conc = {0.1, 0.1, 0.1, 0.1};
MapIndexed[(data[[#2, All, 2]] /= #1) &, conc];

Now solve the linear system. Here I'm using the definition of the linear lease squares just for fun.

soln = With[{A = data[[1 ;; 4, All, 2]]}, 
  Solve[A.Transpose@A.{c1, c2, c3, c4} == A.data[[5, All, 2]], 
  {c1, c2, c3, c4}]]
(* {{c1 -> 0.059183, c2 -> 0.00586333, c3 -> 0.0171604, c4 -> 0.0195068}} *)

(* Note, I would have liked to replace {c1, c2, c3, c4} with the actual filenames but they are not valid mathematica symbols *)

The above expression can also (should also?) be entered as LeastSquares[Transpose@data[[1 ;; 4, All, 2]], data[[5, All, 2]]] which yields the same result as a list of coefficients rather than a list of rules.

Create the theoretical curve and compare it to the original

theoretical = Transpose[{data[[1, All, 1]], 
   {c1, c2, c3, c4}.data[[1 ;; 4, All, 2]] /. First@soln}];
ListLinePlot[{theoretical, data[[5]]}]

Mathematica graphics

If you are interested in parameter errors, then use LinearModelFit.

lm = LinearModelFit[{Transpose@data[[1 ;; 4, All, 2]], data[[5, All, 2]]}];
lm["ParameterTable"]

Mathematica graphics

Not surprisingly, the third variable, corresponding to PT, has the largest error.

share|improve this answer
    
you are right about the PT. The csv file is still saved on the instrument. I couldn't get into the room to retrieve the spectrum. Thank you so much. I'll repost tomorrow with the PT included so that we can get a good answer. –  olliepower Feb 23 at 3:17
    
I consulted my notebook. The PT solution was clear at at this PH. Could you elaborate on the theoretical and actual fit plot? I'm not sure I understand what is being plotted. Thank you so much. I'll repost tomorrow with all of the data processed in case you are interested. –  olliepower Feb 23 at 3:49
    
@olliepower the theoretical plot is simulated from the molar absorptivities of the standard compounds and the best fit concentrations of these standards that makes up the final solution. The other line is the experimental plot. –  bobthechemist Feb 23 at 14:30

Here is the final result if you are interested.

 data = Import /@ 
   FileNames["/home/marco/LatexDocs/analytical/uvvis/PHOS_*.CSV"];
data = data[[All, 2 ;;, {1, 2}]];
ListLinePlot[data, PlotLegends -> {"BB", "MR", "PT", "TB", "UI"}, 
 AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16], 
   Style["Absorbance", FontSize -> 16]}, 
 PlotRange -> {{400, 750}, {0, 0.9}}]
concBB = 1.6*10^(-5);
concMR = 37.1  * 10^(-5);
concPT = 31.4 * 10^(-5);
concTB =  21.4 * 10^(-5);
conc = {concBB, concMR, concPT, concTB};
MapIndexed[(data[[#2, All, 2]] /= #1) &, conc];
soln = With[{A = data[[1 ;; 4, All, 2]]}, 
   Solve[A.Transpose@A.{c1, c2, c3, c4} == A.data[[5, All, 2]], {c1, 
     c2, c3, c4}]];

theoretical = 
  Transpose[{data[[1, All, 
      1]], {c1, c2, c3, c4}.data[[1 ;; 4, All, 2]] /. First@soln}];
ListLinePlot[{theoretical, data[[5]]}, 
 AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16], 
   Style["Absorbance", FontSize -> 16]}]

lm = LinearModelFit[{Transpose@data[[1 ;; 4, All, 2]], 
    data[[5, All, 2]]}];
error1 = Total[(#^2) & @  lm["FitResiduals"]];




data2 = Import /@ 
   FileNames["/home/marco/LatexDocs/analytical/uvvis/CARB_*.CSV"];
data2 = data2[[All, 2 ;;, {1, 2}]];
ListLinePlot[data2, PlotLegends -> {"BB", "MR", "PT", "TB", "UI"}, 
 AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16], 
   Style["Absorbance", FontSize -> 16]}, 
 PlotRange -> {{400, 750}, {0, 0.9}}]

MapIndexed[(data2[[#2, All, 2]] /= #1) &, conc];
soln2 = With[{A = data2[[1 ;; 4, All, 2]]}, 
   Solve[A.Transpose@A.{c1, c2, c3, c4} == A.data2[[5, All, 2]], {c1, 
     c2, c3, c4}]];

theoretical = 
  Transpose[{data2[[1, All, 
      1]], {c1, c2, c3, c4}.data2[[1 ;; 4, All, 2]] /. First@soln}];
ListLinePlot[{theoretical, data2[[5]]}]
lm2 = LinearModelFit[{Transpose@data2[[1 ;; 4, All, 2]], 
    data2[[5, All, 2]]}];
error2 = Total[(#^2) & @  lm2["FitResiduals"]];

data3 = Import /@ 
   FileNames["/home/marco/LatexDocs/analytical/uvvis/KHP_*.CSV"];
data3 = data3[[All, 2 ;;, {1, 2}]];
ListLinePlot[data3, PlotLegends -> {"BB", "MR", "PT", "TB", "UI"}, 
 AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16], 
   Style["Absorbance", FontSize -> 16]}, 
 PlotRange -> {{400, 750}, {0, 0.9}}]

MapIndexed[(data3[[#2, All, 2]] /= #1) &, conc];
soln3 = With[{A = data3[[1 ;; 4, All, 2]]}, 
   Solve[A.Transpose@A.{c1, c2, c3, c4} == A.data3[[5, All, 2]], {c1, 
     c2, c3, c4}]];

theoretical = 
  Transpose[{data3[[1, All, 
      1]], {c1, c2, c3, c4}.data3[[1 ;; 4, All, 2]] /. First@soln}];
ListLinePlot[{theoretical, data3[[5]]}, 
 AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16], 
   Style["Absorbance", FontSize -> 16]}]
lm3 = LinearModelFit[{Transpose@data3[[1 ;; 4, All, 2]], 
    data3[[5, All, 2]]}];
error3 = Total[(#^2) & @  lm3["FitResiduals"]];

dataSummary = {c1, c2, c3, c4} /. {soln, soln2, soln3}
errorSummary = {error1, error2, error3}

Mean Concentration Molar

{{0.000469681, 0.00678397, 0.0169567, 0.0012554}}

Plots: each set of three corresponds to: (spectrum, theoretical, plot of residuals) In Phosphate Buffer Theoretical Comparison Phosphate Plot of Residuals Phosphate

In Carbonate Buffer Carbonate Theoretical Carbonate Error

In KHP buffer KHP theoretical KHP error

share|improve this answer
    
@bobthechemist I have updated the final result if you are interested. Is there any reason why the theoretical plot looks so good for the first plot and so bad for the rest of them? Its pretty odd that there is an apparent structure to the residuals as well. –  olliepower Feb 23 at 18:31

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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