1
$\begingroup$

I am trying to write my own Allan Deviation calculator in Mathematica using the definition, $$\sigma(\tau)^{2}=\frac{1}{2M-2}\sum_{i=1}^{N-1}(\bar{y}_{i+1} - \bar{y}_{i})^{2}\text{.}$$ Where the $\bar{y}$ are mean partitioned data for some multiple of measurement time $\tau$, and $M$ is the length of the number of $(\bar{y}_{i+1} - \bar{y}_{i})^{2}$ associated with each $\tau$. I made a toy model of partitioning the data that I am happy with, but I am struggling to extending it to the next step, which is simply summing the $(\bar{y}_{i+1} - \bar{y}_{i})^{2}$ for each $\tau$. The toy model is:

NumberPoints = 64;
IndexData = Range [NumberPoints];

Do[
  Print["________________τ=", j];
  Do[
   Print[
    "(MEAN", Partition[IndexData, j][[2 i]], "-", "MEAN", 
    Partition[IndexData, j][[2 i - 1]], 
    "\!\(\*SuperscriptBox[\()\), \(2\)]\)"
    ],
   {i, 1, (NumberPoints/(j*2))}
   ],
  {j, 1, NumberPoints/4}
  ];

If I try, for the next step:

Do[
  Clear[BlockMeanDifferenceList];
  BlockMeanDifferenceList = {};
  Do[
   BlockMeanDifferenceList = (
     Mean[Partition[RandomData, τ][[2 i]]] - 
     Mean[Partition[RandomData, τ][[2 i - 1]]]
   )^2;
   AppendTo[σ, (1/(2 Length[BlockMeanDifference] - 2)*Total[BlockMeanDifference])^(1/2)];,
   {i, 1, (NumberOfPoints/(τ*2))}
   ];,
  {τ, 1, NumberOfPoints/4}
  ];

My sigma list is returned for zero in all elements...

Any thoughts?

$\endgroup$
  • $\begingroup$ Your code seems to be missing a few things, e.g. the initialization of σ. Also, what's the difference between BlockMeanDifferenceList and BlockMeanDifference, and NumberOfPoints vs NumberPoints (I assume they're the same?)? Instead of Clear inside a loop, consider using Module. And finally, your code seems to contradict your formula: One has differences of the form $y_2-y_1,y_3-y_2,...$ while the other one has $y_2-y_1,y_4-y_3$ (i.e. disjoint) $\endgroup$ – Lukas Lang Sep 9 '17 at 15:23
  • $\begingroup$ @Mathe172 Sorry, you are right, σ={} should be included. and NumberOfPoints and NumberPoints is the same thing -- I was just sloppy. As for the disjoint, I don't see where that is as the indexing is the same for both the toy model showing the partitioning and for the real attempt? $\endgroup$ – Q.P. Sep 9 '17 at 15:30
  • $\begingroup$ I meant the difference between code and formula (the one you provide), not between the two code snippets $\endgroup$ – Lukas Lang Sep 9 '17 at 15:31
1
$\begingroup$

Below an attempt to implement your formula as I understand it (the prefactor is missing, as I'm not quite sure what $M$ should be) - please do correct me if I misunderstood something.

avar[data_, τ_] := Total[
  (
    Subtract @@@ (*get the difference of the pairs*)
     Partition[ (*create overlapping subsequences of length 2*)
      Mean /@ Partition[data, τ], (*partition data into block and calculate mean*)
      2, 1
      ]
    )^2
  ]

With this, avar[Array[Subscript[a, #] &, 16], 4] yields:

Output of function

If you don't want overlapping differences, simply remove the ,1 in the outer partition (this controls the offset of the individual sublists)

Update

Since you tagged this , here is an approach more closely based on yours (it's impossible to tell you what's wrong with your attempt specifically, as there are still a few obvious typos - see my comment)

avarProc[data_, τ_] := Module[
  {
   partitions = Partition[data, τ],
   σ = {},
   diff
   },
  Do[
   diff = (Mean@partitions[[i + 1]] - Mean@partitions[[i]])^2;
   AppendTo[σ, diff],
   {i, 1, Length@data/τ - 1}
   ];
  Total@σ
  ]

This produces the same output as the above. Again, this uses overlapping sublists. To change this, simply take your indexing and loop bounds.

$\endgroup$
  • $\begingroup$ Thanks for the answer! Can you explain a little what Module does? I have not come across this yet $\endgroup$ – Q.P. Sep 9 '17 at 15:57
  • 1
    $\begingroup$ It basically allows you to introduce local variables. For a more detailed (and accurate) description, see this great answer and the documenation here (Module) and here (Scoping constructs) $\endgroup$ – Lukas Lang Sep 9 '17 at 16:02

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.